From 12e78b79d3d5644e81621103b7411387fd91821b Mon Sep 17 00:00:00 2001 From: Markus Frank <Markus.Frank@cern.ch> Date: Sun, 24 Jul 2022 03:05:26 +0200 Subject: [PATCH] Support replicates and parameterised volumes when converting to Geant4 --- DDCore/include/DD4hep/BitField64.h | 2 +- DDCore/include/DD4hep/BitFieldCoder.h | 2 +- DDCore/include/DD4hep/CartesianGridXY.h | 2 +- DDCore/include/DD4hep/ComponentProperties.h | 2 +- DDCore/include/DD4hep/DD4hepRootPersistency.h | 2 +- DDCore/include/DD4hep/Detector.h | 28 +- DDCore/include/DD4hep/Factories.h | 6 +- DDCore/include/DD4hep/IDDescriptor.h | 7 +- DDCore/include/DD4hep/VolumeManager.h | 12 +- DDCore/include/DD4hep/Volumes.h | 258 +++++++++++++-- DDCore/include/Parsers/detail/Dimension.h | 4 + DDCore/include/Parsers/detail/Dimension.imp | 2 + DDCore/include/XML/UnicodeValues.h | 1 + DDCore/src/GeoDictionary.h | 2 + DDCore/src/IDDescriptor.cpp | 13 +- DDCore/src/Volumes.cpp | 313 +++++++++++++++--- DDG4/include/DDG4/Geant4GeometryInfo.h | 22 +- DDG4/include/DDG4/Geant4Helpers.h | 88 +++++ .../DDG4/Geant4PlacementParameterisation.h | 98 ++++++ DDG4/include/DDG4/Geant4SensDetAction.h | 37 +-- DDG4/include/DDG4/Geant4SensDetAction.inl | 26 +- DDG4/include/DDG4/Geant4Vertex.h | 2 +- DDG4/include/DDG4/Geant4VolumeManager.h | 17 +- DDG4/src/Geant4Converter.cpp | 122 ++++--- DDG4/src/Geant4Helpers.cpp | 199 +++++++++++ DDG4/src/Geant4InputHandling.cpp | 22 +- DDG4/src/Geant4Particle.cpp | 22 +- DDG4/src/Geant4ParticleHandler.cpp | 38 +-- DDG4/src/Geant4PlacementParameterisation.cpp | 102 ++++++ DDG4/src/Geant4SensDetAction.cpp | 38 +-- DDG4/src/Geant4ShapeConverter.cpp | 88 +++-- DDG4/src/Geant4VolumeManager.cpp | 144 ++++++-- DDG4/tpython/DDPython.cpp | 4 + .../ClientTests/compact/ParamVolume1D.xml | 162 +++++++++ .../ClientTests/compact/ParamVolume2D.xml | 83 +++++ .../ClientTests/compact/ParamVolume3D.xml | 95 ++++++ examples/ClientTests/scripts/ParamVolume.mac | 11 + examples/ClientTests/scripts/ParamVolume.py | 124 +++++++ .../ClientTests/scripts/SiliconBlockGFlash.py | 36 +- examples/ClientTests/src/ParamVolume_geo.cpp | 99 ++++++ 40 files changed, 1989 insertions(+), 346 deletions(-) create mode 100644 DDG4/include/DDG4/Geant4Helpers.h create mode 100644 DDG4/include/DDG4/Geant4PlacementParameterisation.h create mode 100644 DDG4/src/Geant4Helpers.cpp create mode 100644 DDG4/src/Geant4PlacementParameterisation.cpp create mode 100644 examples/ClientTests/compact/ParamVolume1D.xml create mode 100644 examples/ClientTests/compact/ParamVolume2D.xml create mode 100644 examples/ClientTests/compact/ParamVolume3D.xml create mode 100644 examples/ClientTests/scripts/ParamVolume.mac create mode 100644 examples/ClientTests/scripts/ParamVolume.py create mode 100644 examples/ClientTests/src/ParamVolume_geo.cpp diff --git a/DDCore/include/DD4hep/BitField64.h b/DDCore/include/DD4hep/BitField64.h index 2b77f34cf..5e4af0288 100644 --- a/DDCore/include/DD4hep/BitField64.h +++ b/DDCore/include/DD4hep/BitField64.h @@ -14,7 +14,7 @@ #define DD4HEP_BITFIELD64_H // Framework include files -#include "DDSegmentation/BitField64.h" +#include <DDSegmentation/BitField64.h> /// Namespace for the AIDA detector description toolkit namespace dd4hep { diff --git a/DDCore/include/DD4hep/BitFieldCoder.h b/DDCore/include/DD4hep/BitFieldCoder.h index a628cfdfd..96a985b90 100644 --- a/DDCore/include/DD4hep/BitFieldCoder.h +++ b/DDCore/include/DD4hep/BitFieldCoder.h @@ -14,7 +14,7 @@ #define DD4HEP_BITFIELDCODER_H // Framework include files -#include "DDSegmentation/BitFieldCoder.h" +#include <DDSegmentation/BitFieldCoder.h> /// Namespace for the AIDA detector description toolkit namespace dd4hep { diff --git a/DDCore/include/DD4hep/CartesianGridXY.h b/DDCore/include/DD4hep/CartesianGridXY.h index 0107d0ae3..2b8057a5e 100644 --- a/DDCore/include/DD4hep/CartesianGridXY.h +++ b/DDCore/include/DD4hep/CartesianGridXY.h @@ -16,7 +16,7 @@ #define DD4HEP_CARTESIANGRIDXY_H 1 // Framework include files -#include "DD4hep/Segmentations.h" +#include <DD4hep/Segmentations.h> /// Namespace for the AIDA detector description toolkit namespace dd4hep { diff --git a/DDCore/include/DD4hep/ComponentProperties.h b/DDCore/include/DD4hep/ComponentProperties.h index aefcbf4e7..bb354cb40 100644 --- a/DDCore/include/DD4hep/ComponentProperties.h +++ b/DDCore/include/DD4hep/ComponentProperties.h @@ -14,7 +14,7 @@ #define DD4HEP_COMPONENTPROPERTIES_H // Framework include files -#include "DD4hep/Grammar.h" +#include <DD4hep/Grammar.h> // C/C++ include files #include <algorithm> diff --git a/DDCore/include/DD4hep/DD4hepRootPersistency.h b/DDCore/include/DD4hep/DD4hepRootPersistency.h index 45d047a02..edf04cd34 100644 --- a/DDCore/include/DD4hep/DD4hepRootPersistency.h +++ b/DDCore/include/DD4hep/DD4hepRootPersistency.h @@ -14,7 +14,7 @@ #define DD4HEP_DD4HEPROOTPERSISTENCY_H // Framework include files -#include "DD4hep/DetectorData.h" +#include <DD4hep/DetectorData.h> /// Helper class to support ROOT persistency of Detector objects diff --git a/DDCore/include/DD4hep/Detector.h b/DDCore/include/DD4hep/Detector.h index d9b988b28..46a4c0d95 100644 --- a/DDCore/include/DD4hep/Detector.h +++ b/DDCore/include/DD4hep/Detector.h @@ -13,22 +13,22 @@ #ifndef DD4HEP_DETECTOR_H #define DD4HEP_DETECTOR_H -#include "DD4hep/Version.h" +#include <DD4hep/Version.h> // Framework includes -#include "DD4hep/Handle.h" -#include "DD4hep/Fields.h" -#include "DD4hep/Objects.h" -#include "DD4hep/Shapes.h" -#include "DD4hep/Volumes.h" -#include "DD4hep/Readout.h" -#include "DD4hep/DetElement.h" -#include "DD4hep/NamedObject.h" -#include "DD4hep/Segmentations.h" -#include "DD4hep/VolumeManager.h" -#include "DD4hep/OpticalSurfaceManager.h" -#include "DD4hep/ExtensionEntry.h" -#include "DD4hep/BuildType.h" +#include <DD4hep/Handle.h> +#include <DD4hep/Fields.h> +#include <DD4hep/Objects.h> +#include <DD4hep/Shapes.h> +#include <DD4hep/Volumes.h> +#include <DD4hep/Readout.h> +#include <DD4hep/DetElement.h> +#include <DD4hep/NamedObject.h> +#include <DD4hep/Segmentations.h> +#include <DD4hep/VolumeManager.h> +#include <DD4hep/OpticalSurfaceManager.h> +#include <DD4hep/ExtensionEntry.h> +#include <DD4hep/BuildType.h> // C/C++ include files #include <map> diff --git a/DDCore/include/DD4hep/Factories.h b/DDCore/include/DD4hep/Factories.h index 1fb780f15..8caf059b0 100644 --- a/DDCore/include/DD4hep/Factories.h +++ b/DDCore/include/DD4hep/Factories.h @@ -14,9 +14,9 @@ #define DD4HEP_FACTORIES_H // Framework include files -#include "DD4hep/Plugins.h" -#include "DD4hep/DetElement.h" -#include "DD4hep/NamedObject.h" +#include <DD4hep/Plugins.h> +#include <DD4hep/DetElement.h> +#include <DD4hep/NamedObject.h> // C/C++ include files #include <cstdarg> diff --git a/DDCore/include/DD4hep/IDDescriptor.h b/DDCore/include/DD4hep/IDDescriptor.h index f691196ec..2e33c80c1 100644 --- a/DDCore/include/DD4hep/IDDescriptor.h +++ b/DDCore/include/DD4hep/IDDescriptor.h @@ -36,8 +36,9 @@ namespace dd4hep { */ class IDDescriptor: public Handle<IDDescriptorObject> { public: - typedef std::vector<std::pair<std::string, const BitFieldElement*> > FieldMap; - typedef std::vector<std::pair<size_t, std::string> > FieldIDs; + typedef BitFieldElement Field; + typedef std::vector<std::pair<std::string, const Field*> > FieldMap; + typedef std::vector<std::pair<size_t, std::string> > FieldIDs; public: /// Default constructor @@ -61,6 +62,8 @@ namespace dd4hep { /// Get the field descriptor of one field by its identifier const BitFieldElement* field(size_t identifier) const; #ifndef __MAKECINT__ + /// Encode partial volume identifiers to a volumeID. + static VolumeID encode(const Field* fld, VolumeID value); /// Encode a set of volume identifiers (corresponding to this description of course!) to a volumeID. VolumeID encode(const std::vector<std::pair<std::string, int> >& ids) const; /// Encode a set of volume identifiers to a volumeID with the system ID on the top bits diff --git a/DDCore/include/DD4hep/VolumeManager.h b/DDCore/include/DD4hep/VolumeManager.h index 0e21d5160..9a477813b 100644 --- a/DDCore/include/DD4hep/VolumeManager.h +++ b/DDCore/include/DD4hep/VolumeManager.h @@ -15,14 +15,14 @@ #define DD4HEP_VOLUMEMANAGER_H // Framework include files -#include "DD4hep/Volumes.h" -#include "DD4hep/DetElement.h" -#include "DD4hep/NamedObject.h" -#include "DD4hep/IDDescriptor.h" -#include "DD4hep/ConditionsMap.h" +#include <DD4hep/Volumes.h> +#include <DD4hep/DetElement.h> +#include <DD4hep/NamedObject.h> +#include <DD4hep/IDDescriptor.h> +#include <DD4hep/ConditionsMap.h> // ROOT include files -#include "TGeoMatrix.h" +#include <TGeoMatrix.h> /// Namespace for the AIDA detector description toolkit namespace dd4hep { diff --git a/DDCore/include/DD4hep/Volumes.h b/DDCore/include/DD4hep/Volumes.h index cab011a43..98ea1231e 100644 --- a/DDCore/include/DD4hep/Volumes.h +++ b/DDCore/include/DD4hep/Volumes.h @@ -14,20 +14,19 @@ #define DD4HEP_VOLUMES_H // Framework include files -#include "DD4hep/Handle.h" -#include "DD4hep/Shapes.h" -#include "DD4hep/Objects.h" +#include <DD4hep/Handle.h> +#include <DD4hep/Shapes.h> +#include <DD4hep/Objects.h> +#include <DD4hep/BitField64.h> // C/C++ include files #include <map> #include <memory> // ROOT include file (includes TGeoVolume + TGeoShape) -#include "TGeoNode.h" -#include "TGeoPatternFinder.h" - -// Recent ROOT versions -#include "TGeoExtension.h" +#include <TGeoNode.h> +#include <TGeoPatternFinder.h> +#include <TGeoExtension.h> /// Namespace for the AIDA detector description toolkit namespace dd4hep { @@ -116,12 +115,64 @@ namespace dd4hep { /// String representation for debugging std::string str() const; }; + /// Optional parameters to implement special features such as parametrization + /** + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_CORE + */ + class Parameterisation { + public: + /** Feature: Geant4 parameterised volumes */ + using Dimension = std::pair<Transform3D, size_t>; + /// Reference to the starting point of the parameterisation + Transform3D start { }; + /// Reference to the parameterised transformation for dimension 1 + Dimension trafo1D { {}, 0UL }; + /// Reference to the parameterised transformation for dimension 2 + Dimension trafo2D { {}, 0UL }; + /// Reference to the parameterised transformation for dimension 3 + Dimension trafo3D { {}, 0UL }; + /// Number of entries for the parameterisation in dimension 2 + unsigned long flags { 0 }; + /// Number of entries for the parameterisation in dimension 2 + unsigned long refCount { 0 }; + /// Reference to the placements of this volume + std::vector<PlacedVolume> placements { }; + /// Bitfield from sensitive detector to encode the volume ID on the fly + const detail::BitFieldElement* field { nullptr }; + + public: + /// Default constructor + Parameterisation() = default; + /// Default move + Parameterisation(Parameterisation&& copy) = default; + /// Copy constructor + Parameterisation(const Parameterisation& c) = default; + /// Default destructor + virtual ~Parameterisation() = default; + /// Move assignment + Parameterisation& operator=(Parameterisation&& copy) = default; + /// Assignment operator + Parameterisation& operator=(const Parameterisation& copy) = default; + /// Increase ref count + Parameterisation* addref(); + /// Decrease ref count + void release(); + /// Enable ROOT persistency + ClassDef(Parameterisation,200); + }; /// Magic word to detect memory corruptions - unsigned long magic = 0; + unsigned long magic { 0 }; /// Reference count on object (used to implement Grab/Release) - long refCount = 0; + long refCount { 0 }; + /// Reference to the parameterised transformation + Parameterisation* params { nullptr }; /// ID container VolIDs volIDs; + + public: /// Default constructor PlacedVolumeExtension(); /// Default move @@ -133,13 +184,15 @@ namespace dd4hep { /// Move assignment PlacedVolumeExtension& operator=(PlacedVolumeExtension&& copy) { magic = std::move(copy.magic); + params = std::move(copy.params); volIDs = std::move(copy.volIDs); return *this; } /// Assignment operator - PlacedVolumeExtension& operator=(const PlacedVolumeExtension& c) { - magic = c.magic; - volIDs = c.volIDs; + PlacedVolumeExtension& operator=(const PlacedVolumeExtension& copy) { + magic = copy.magic; + params = copy.params; + volIDs = copy.volIDs; return *this; } /// TGeoExtension overload: Method called whenever requiring a pointer to the extension @@ -149,6 +202,10 @@ namespace dd4hep { /// Enable ROOT persistency ClassDefOverride(PlacedVolumeExtension,200); }; + class PlacedVolumeFeatureExtension : public PlacedVolumeExtension { + public: + public: + }; /// Handle class holding a placed volume (also called physical volume) /** @@ -201,6 +258,8 @@ namespace dd4hep { Object* data() const; /// Access the copy number of this placement within its mother int copyNumber() const; + /// Access placed volume flags + long flags() const; /// Volume material Material material() const; /// Logical volume of this placement @@ -311,8 +370,18 @@ namespace dd4hep { VETO_SIMU = 1, VETO_RECO = 2, VETO_DISPLAY = 3, - REFLECTED = 10 + REFLECTED = 10, + }; + enum ReplicationAxis { + REPLICATED = 1UL << 4, + PARAMETERIZED = 1UL << 5, + X_axis = 1UL << 8, + Y_axis = 1UL << 9, + Z_axis = 1UL << 10, + Rho_axis = 1UL << 11, + Phi_axis = 1UL << 12 }; + public: /// Default constructor Volume() = default; @@ -388,35 +457,162 @@ namespace dd4hep { PlacedVolume placeVolume(const Volume& volume, TGeoMatrix* tr) const; /// Place daughter volume with generic TGeo matrix PlacedVolume placeVolume(const Volume& volume, int copy_nr, TGeoMatrix* tr) const; - /// Parametrized volume implementation - /** Embedding parametrized daughter placements in a mother volume - * @param start start transormation for the first placement - * @param count Number of entities to be placed + /// 1D volume replication implementation + /** Embedding parameterised daughter placements in a mother volume * @param entity Daughter volume to be placed + * @param axis Replication axis direction in the frame of the mother + * @param count Number of entities to be placed * @param inc Transformation increment for each iteration + * @param start start transormation for the first placement */ - void paramVolume1D(const Transform3D& start, size_t count, Volume entity, const Transform3D& inc); - /// Parametrized volume implementation - /** Embedding parametrized daughter placements in a mother volume - * @param count Number of entities to be placed + PlacedVolume replicate(const Volume entity, ReplicationAxis axis, size_t count, double inc, double start=0e0); + /// 1D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param start start transormation for the first placement * @param entity Daughter volume to be placed + * @param count Number of entities to be placed * @param inc Transformation increment for each iteration */ - void paramVolume1D(size_t count, Volume entity, const Transform3D& trafo); - /// Parametrized volume implementation - /** Embedding parametrized daughter placements in a mother volume - * @param count Number of entities to be placed + PlacedVolume paramVolume1D(const Transform3D& start, Volume entity, size_t count, const Transform3D& inc); + /// 1D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume * @param entity Daughter volume to be placed + * @param count Number of entities to be placed * @param inc Transformation increment for each iteration */ - void paramVolume1D(size_t count, Volume entity, const Position& inc); - /// Parametrized volume implementation - /** Embedding parametrized daughter placements in a mother volume + PlacedVolume paramVolume1D(Volume entity, size_t count, const Transform3D& trafo); + /// 1D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param entity Daughter volume to be placed * @param count Number of entities to be placed + * @param inc Transformation increment for each iteration + */ + PlacedVolume paramVolume1D(Volume entity, size_t count, const Position& inc); + /// 1D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume * @param entity Daughter volume to be placed + * @param count Number of entities to be placed * @param inc Transformation increment for each iteration */ - void paramVolume1D(size_t count, Volume entity, const RotationZYX& inc); + PlacedVolume paramVolume1D(Volume entity, size_t count, const RotationZYX& inc); + + /// 2D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + */ + PlacedVolume paramVolume2D(Volume entity, + size_t count_1, const Transform3D& inc_1, + size_t count_2, const Transform3D& inc_2); + + /// Constructor to be used when creating a new parameterised volume object + /** Embedding parameterised daughter placements in a mother volume + * @param start start transormation for the first placement + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + */ + PlacedVolume paramVolume2D(const Transform3D& start, + Volume entity, + size_t count_1, + const Position& inc_1, + size_t count_2, + const Position& inc_2); + + /// Constructor to be used when creating a new parameterised volume object + /** Embedding parameterised daughter placements in a mother volume + * @param start start transormation for the first placement + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + */ + PlacedVolume paramVolume2D(Volume entity, + size_t count_1, + const Position& inc_1, + size_t count_2, + const Position& inc_2); + + /// 2D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param start start transormation for the first placement + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + */ + PlacedVolume paramVolume2D(const Transform3D& start, Volume entity, + size_t count_1, const Transform3D& inc_1, + size_t count_2, const Transform3D& inc_2); + + /// 3D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + * @param count_3 Number of entities to be placed in dimension 3 + * @param inc_3 Transformation increment for each iteration in dimension 3 + */ + PlacedVolume paramVolume3D(Volume entity, + size_t count_1, const Transform3D& inc_1, + size_t count_2, const Transform3D& inc_2, + size_t count_3, const Transform3D& inc_3); + + /// 3D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param start start transormation for the first placement + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + * @param count_3 Number of entities to be placed in dimension 3 + * @param inc_3 Transformation increment for each iteration in dimension 3 + */ + PlacedVolume paramVolume3D(const Transform3D& start, Volume entity, + size_t count_1, const Transform3D& inc_1, + size_t count_2, const Transform3D& inc_2, + size_t count_3, const Transform3D& inc_3); + + /// 3D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + * @param count_3 Number of entities to be placed in dimension 3 + * @param inc_3 Transformation increment for each iteration in dimension 3 + */ + PlacedVolume paramVolume3D(Volume entity, + size_t count_1, const Position& inc_1, + size_t count_2, const Position& inc_2, + size_t count_3, const Position& inc_3); + + /// 3D Parameterised volume implementation + /** Embedding parameterised daughter placements in a mother volume + * @param start start transormation for the first placement + * @param entity Daughter volume to be placed + * @param count_1 Number of entities to be placed in dimension 1 + * @param inc_1 Transformation increment for each iteration in dimension 1 + * @param count_2 Number of entities to be placed in dimension 2 + * @param inc_2 Transformation increment for each iteration in dimension 2 + * @param count_3 Number of entities to be placed in dimension 3 + * @param inc_3 Transformation increment for each iteration in dimension 3 + */ + PlacedVolume paramVolume3D(const Transform3D& start, Volume entity, + size_t count_1, const Position& inc_1, + size_t count_2, const Position& inc_2, + size_t count_3, const Position& inc_3); /// Set user flags in bit-field void setFlagBit(unsigned int bit); @@ -427,7 +623,7 @@ namespace dd4hep { bool isReflected() const; /// Set the volume's option value - void setOption(const std::string& opt) const; + const Volume& setOption(const std::string& opt) const; /// Access the volume's option value std::string option() const; diff --git a/DDCore/include/Parsers/detail/Dimension.h b/DDCore/include/Parsers/detail/Dimension.h index 93fd22715..98b4ad2d7 100644 --- a/DDCore/include/Parsers/detail/Dimension.h +++ b/DDCore/include/Parsers/detail/Dimension.h @@ -685,6 +685,10 @@ namespace dd4hep { Dimension position(bool throw_if_not_present = true) const; /// Access child element with tag "rotation" as Dimension object Dimension rotation(bool throw_if_not_present = true) const; + /// Access child element with tag "transformation" as Dimension object + Dimension transformation(bool throw_if_not_present = true) const; + /// Access child element with tag "transform" as Dimension object + Dimension transform(bool throw_if_not_present = true) const; /// Access child element with tag "cone" as Dimension object Dimension cone(bool throw_if_not_present = true) const; /// Access child element with tag "sphere" as Dimension object diff --git a/DDCore/include/Parsers/detail/Dimension.imp b/DDCore/include/Parsers/detail/Dimension.imp index d7dff15e5..f22e6cc03 100644 --- a/DDCore/include/Parsers/detail/Dimension.imp +++ b/DDCore/include/Parsers/detail/Dimension.imp @@ -222,6 +222,8 @@ XML_CHILD_ACCESSOR_XML_DIM(params) XML_CHILD_ACCESSOR_XML_DIM(parameters) XML_CHILD_ACCESSOR_XML_DIM(position) XML_CHILD_ACCESSOR_XML_DIM(rotation) +XML_CHILD_ACCESSOR_XML_DIM(transform) +XML_CHILD_ACCESSOR_XML_DIM(transformation) XML_CHILD_ACCESSOR_XML_DIM(cone) XML_CHILD_ACCESSOR_XML_DIM(sphere) XML_CHILD_ACCESSOR_XML_DIM(torus) diff --git a/DDCore/include/XML/UnicodeValues.h b/DDCore/include/XML/UnicodeValues.h index 576ff45dd..0e088210f 100644 --- a/DDCore/include/XML/UnicodeValues.h +++ b/DDCore/include/XML/UnicodeValues.h @@ -497,6 +497,7 @@ UNICODE (torus); UNICODE (tracker); UNICODE (tracking_cylinder); UNICODE (tracking_volume); +UNICODE (transform); UNICODE (transformation); UNICODE (trap); UNICODE (trapezoid); diff --git a/DDCore/src/GeoDictionary.h b/DDCore/src/GeoDictionary.h index d1f459bd9..52a3ef11a 100644 --- a/DDCore/src/GeoDictionary.h +++ b/DDCore/src/GeoDictionary.h @@ -55,6 +55,8 @@ template vector<pair<string, int> >::iterator; #pragma link C++ class vector<pair<string, int> >+; #pragma link C++ class vector<pair<string, int> >::iterator; #pragma link C++ class dd4hep::PlacedVolumeExtension::VolIDs+; +#pragma link C++ class dd4hep::PlacedVolumeExtension::Parameterisation+; +#pragma link C++ class dd4hep::PlacedVolumeExtension::Parameterisation::Dimension+; #pragma link C++ class dd4hep::PlacedVolumeExtension+; #pragma link C++ class vector<dd4hep::PlacedVolume>+; #pragma link C++ class dd4hep::Handle<TGeoNode>+; diff --git a/DDCore/src/IDDescriptor.cpp b/DDCore/src/IDDescriptor.cpp index 7014fd44b..4c163d139 100644 --- a/DDCore/src/IDDescriptor.cpp +++ b/DDCore/src/IDDescriptor.cpp @@ -127,8 +127,7 @@ VolumeID IDDescriptor::get_mask(const std::vector<std::pair<std::string, int> >& } /// Encode a set of volume identifiers (corresponding to this description of course!) to a volumeID. -VolumeID IDDescriptor::encode(const std::vector<std::pair<std::string, int> >& id_vector) const -{ +VolumeID IDDescriptor::encode(const std::vector<std::pair<std::string, int> >& id_vector) const { VolumeID id = 0; //const PlacedVolume::VolIDs* ids = (const PlacedVolume::VolIDs*)&id_vector; //printout(INFO,"IDDescriptor","VolIDs: %s",ids->str().c_str()); @@ -141,6 +140,16 @@ VolumeID IDDescriptor::encode(const std::vector<std::pair<std::string, int> >& i return id; } +/// Encode partial volume identifiers to a volumeID. +VolumeID IDDescriptor::encode(const Field* fld, VolumeID value) { + if ( fld ) { + int off = fld->offset(); + return ((fld->value(value << off) << off)&fld->mask()); + } + except("IDDescriptor","dd4hep: %s: Cannot encode value with void Field reference."); + return 0UL; +} + /// Encode a set of volume identifiers to a volumeID with the system ID on the top bits VolumeID IDDescriptor::encode_reverse(const std::vector<std::pair<std::string, int> >& id_vector) const { diff --git a/DDCore/src/Volumes.cpp b/DDCore/src/Volumes.cpp index 415963e68..00a38f6f4 100644 --- a/DDCore/src/Volumes.cpp +++ b/DDCore/src/Volumes.cpp @@ -12,25 +12,26 @@ //========================================================================== // Framework include files -#include "DD4hep/Detector.h" -#include "DD4hep/Printout.h" -#include "DD4hep/InstanceCount.h" -#include "DD4hep/MatrixHelpers.h" -#include "DD4hep/detail/ObjectsInterna.h" +#include <DD4hep/Detector.h> +#include <DD4hep/Printout.h> +#include <DD4hep/InstanceCount.h> +#include <DD4hep/MatrixHelpers.h> +#include <DD4hep/detail/ObjectsInterna.h> // ROOT include files -#include "TClass.h" -#include "TColor.h" -#include "TGeoShape.h" -#include "TGeoVolume.h" -#include "TGeoNode.h" -#include "TGeoMatrix.h" -#include "TGeoMedium.h" - -#include "TGeoVoxelFinder.h" -#include "TGeoShapeAssembly.h" -#include "TGeoScaledShape.h" -#include "TMap.h" +#include <TROOT.h> +#include <TClass.h> +#include <TColor.h> +#include <TGeoShape.h> +#include <TGeoVolume.h> +#include <TGeoNode.h> +#include <TGeoMatrix.h> +#include <TGeoMedium.h> + +#include <TGeoVoxelFinder.h> +#include <TGeoShapeAssembly.h> +#include <TGeoScaledShape.h> +#include <TMap.h> // C/C++ include files #include <climits> @@ -54,7 +55,10 @@ typedef TGeoNode geo_node_t; typedef TGeoVolume geo_volume_t; typedef TGeoVolumeAssembly geo_assembly_t; +/// Enable ROOT persistency +ClassImp(VolumeExtension) ClassImp(PlacedVolumeExtension) +ClassImp(PlacedVolumeExtension::Parameterisation) namespace { @@ -336,6 +340,17 @@ PlacedVolume::Processor::Processor() { PlacedVolume::Processor::~Processor() { } +/// Increase ref count +PlacedVolumeExtension::Parameterisation* PlacedVolumeExtension::Parameterisation::addref() { + ++refCount; + return this; +} +/// Decrease ref count +void PlacedVolumeExtension::Parameterisation::release() { + --refCount; + if ( 0 == refCount ) delete this; +} + /// Default constructor PlacedVolumeExtension::PlacedVolumeExtension() : TGeoExtension(), volIDs() @@ -360,6 +375,7 @@ PlacedVolumeExtension::PlacedVolumeExtension(const PlacedVolumeExtension& c) /// Default destructor PlacedVolumeExtension::~PlacedVolumeExtension() { + if ( this->params ) this->params->release(); DECREMENT_COUNTER; } @@ -447,7 +463,26 @@ const PlacedVolume::VolIDs& PlacedVolume::volIDs() const { /// Add identifier PlacedVolume& PlacedVolume::addPhysVolID(const string& nam, int value) { - _data(*this)->volIDs.emplace_back(nam, value); + auto* o = _data(*this); + if ( !o->params ) { + o->volIDs.emplace_back(nam, value); + return *this; + } + if ( value > 0 ) { + except("PlacedVolume", + "+++ addPhysVolID(%s): parameterised volumes only accept '0' is volID." + "These automatically get overwritten with the copy number!", + ptr()->GetName()); + } + if ( !o->volIDs.empty() ) { + except("PlacedVolume", + "+++ addPhysVolID(%s): parameterised volumes can only host 1 physical volume ID." + " vol id '%s' is already defined!", ptr()->GetName(), o->volIDs[0].first.c_str()); + } + for(PlacedVolume pv : o->params->placements) { + auto* p = _data(pv); + p->volIDs.emplace_back(nam, pv->GetNumber()); + } return *this; } @@ -478,9 +513,6 @@ string PlacedVolume::toString() const { return str.str(); } -/// Enable ROOT persistency -ClassImp(VolumeExtension) - /// Default constructor VolumeExtension::VolumeExtension() : TGeoExtension(), region(), limits(), vis(), sens_det() { @@ -690,7 +722,8 @@ PlacedVolume _addNode(TGeoVolume* par, TGeoVolume* daughter, int id, TGeoMatrix* if ( nam_id != n->GetName() ) { printout(ERROR,"PlacedVolume","++ FAILED to place node %s",(const char*)nam_id); } - n->geo_node_t::SetUserExtension(new PlacedVolume::Object()); + PlacedVolume::Object* extension = new PlacedVolume::Object(); + n->geo_node_t::SetUserExtension(extension); return PlacedVolume(n); } @@ -776,37 +809,231 @@ PlacedVolume Volume::placeVolume(const Volume& volume, int copy_no, const Rotati return _addNode(m_element, volume, copy_no, rot); } -/// Constructor to be used when creating a new parametrized volume object -void Volume::paramVolume1D(size_t count, Volume entity, const Position& inc) { - paramVolume1D(Transform3D(), count, entity, Transform3D(inc)); +/// 1D volume replication implementation +PlacedVolume Volume::replicate(const Volume entity, ReplicationAxis axis, + size_t count, double inc, double start) +{ + Transform3D offset(1e0, 0e0, 0e0, axis == X_axis ? start : 0e0, + 0e0, 1e0, 0e0, axis == Y_axis ? start : 0e0, + 0e0, 0e0, 1e0, axis == Z_axis ? start : 0e0); + Transform3D tr(1e0, 0e0, 0e0, axis == X_axis ? inc : 0e0, + 0e0, 1e0, 0e0, axis == Y_axis ? inc : 0e0, + 0e0, 0e0, 1e0, axis == Z_axis ? inc : 0e0); + PlacedVolume pv = paramVolume1D(offset, entity, count, tr); + auto* data = pv.data(); + data->params->flags = axis | REPLICATED; + return pv; +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume1D(Volume entity, size_t count, const Position& inc) { + return paramVolume1D(Transform3D(), entity, count, Transform3D(inc)); +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume1D(Volume entity, size_t count, const RotationZYX& inc) { + return paramVolume1D(Transform3D(), entity, count, Transform3D(inc)); +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume1D(Volume entity, size_t count, const Transform3D& inc) { + return paramVolume1D(Transform3D(), entity, count, inc); +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume1D(const Transform3D& start, + Volume entity, + size_t count, + const Transform3D& trafo) +{ + Transform3D tr(start); + PlacedVolume pv = + _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr)); + auto* data = pv.data(); + if ( pv->GetNdaughters() > 1 ) { + except("Volume","paramVolume1D: Mother %s has too many daughters: %ld " + "Parameterized volumes may only have one single daughter!", + ptr()->GetName(), pv.volume()->GetNdaughters()); + } + data->params = new PlacedVolumeExtension::Parameterisation(); + data->params->addref(); + data->params->flags = PARAMETERIZED; + data->params->start = start; + data->params->trafo1D.first = trafo; + data->params->trafo1D.second = count; + data->params->trafo2D.second = 0; + data->params->trafo3D.second = 0; + data->params->placements.emplace_back(pv); + for(size_t i=1; i < count; ++i) { + tr *= trafo; + PlacedVolume ppv = + _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr)); + data->params->placements.emplace_back(ppv); + ppv.data()->params = data->params->addref(); + } + return pv; } -/// Constructor to be used when creating a new parametrized volume object -void Volume::paramVolume1D(size_t count, Volume entity, const RotationZYX& inc) { - paramVolume1D(Transform3D(), count, entity, Transform3D(inc)); +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume2D(Volume entity, + size_t count_1, + const Transform3D& trafo_1, + size_t count_2, + const Transform3D& trafo_2) +{ + return paramVolume2D(Transform3D(), entity, count_1, trafo_1, count_2, trafo_2); } -/// Constructor to be used when creating a new parametrized volume object -void Volume::paramVolume1D(size_t count, Volume entity, const Transform3D& inc) { - paramVolume1D(Transform3D(), count, entity, inc); +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume2D(const Transform3D& start, + Volume entity, + size_t count_1, + const Position& pos_1, + size_t count_2, + const Position& pos_2) +{ + return paramVolume2D(start, entity, count_1, Transform3D(pos_1), count_2, Transform3D(pos_2)); } -/// Constructor to be used when creating a new parametrized volume object -void Volume::paramVolume1D(const Transform3D& start, - size_t count, - Volume entity, - const Transform3D& trafo) +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume2D(Volume entity, + size_t count_1, + const Position& pos_1, + size_t count_2, + const Position& pos_2) +{ + return paramVolume2D(Transform3D(), entity, count_1, Transform3D(pos_1), count_2, Transform3D(pos_2)); +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume2D(const Transform3D& start, + Volume entity, + size_t count_1, + const Transform3D& trafo_1, + size_t count_2, + const Transform3D& trafo_2) +{ + PlacedVolume pv = + _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(start)); + auto* data = pv.data(); + if ( pv->GetNdaughters() > 1 ) { + except("Volume","paramVolume1D: Mother %s has too many daughters: %ld " + "Parameterized volumes may only have one single daughter!", + ptr()->GetName(), pv.volume()->GetNdaughters()); + } + data->params = new PlacedVolumeExtension::Parameterisation(); + data->params->addref(); + data->params->flags = PARAMETERIZED; + data->params->start = start; + data->params->trafo1D.first = trafo_1; + data->params->trafo1D.second = count_1; + data->params->trafo2D.first = trafo_2; + data->params->trafo2D.second = count_2; + data->params->trafo3D.second = 0; + data->params->placements.emplace_back(pv); + Transform3D tr2(start); + for(size_t j=0; j < count_2; ++j) { + Transform3D tr1 = tr2; + for(size_t i = 0; i < count_1; ++i) { + if ( !( i == 0 && j == 0 ) ) { + PlacedVolume ppv = + _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr1)); + data->params->placements.emplace_back(ppv); + ppv.data()->params = data->params->addref(); + } + tr1 *= trafo_1; + } + tr2 *= trafo_2; + } + return pv; +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume3D(const Transform3D& start, + Volume entity, + size_t count_1, + const Position& pos_1, + size_t count_2, + const Position& pos_2, + size_t count_3, + const Position& pos_3) { - Transform3D transformation(start); - for(size_t i=0; i<count; ++i) { - _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(transformation)); - transformation *= trafo; + return paramVolume3D(start, entity, + count_1, Transform3D(pos_1), + count_2, Transform3D(pos_2), + count_3, Transform3D(pos_3)); +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume3D(Volume entity, + size_t count_1, + const Position& pos_1, + size_t count_2, + const Position& pos_2, + size_t count_3, + const Position& pos_3) +{ + return paramVolume3D(Transform3D(), entity, + count_1, Transform3D(pos_1), + count_2, Transform3D(pos_2), + count_3, Transform3D(pos_3)); +} + +/// Constructor to be used when creating a new parameterised volume object +PlacedVolume Volume::paramVolume3D(const Transform3D& start, + Volume entity, + size_t count_1, + const Transform3D& trafo_1, + size_t count_2, + const Transform3D& trafo_2, + size_t count_3, + const Transform3D& trafo_3) +{ + PlacedVolume pv = + _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(start)); + auto* data = pv.data(); + if ( pv->GetNdaughters() > 1 ) { + except("Volume","paramVolume1D: Mother %s has too many daughters: %ld " + "Parameterized volumes may only have one single daughter!", + ptr()->GetName(), pv.volume()->GetNdaughters()); } + data->params = new PlacedVolumeExtension::Parameterisation(); + data->params->addref(); + data->params->flags = PARAMETERIZED; + data->params->start = start; + data->params->trafo2D.first = trafo_2; + data->params->trafo2D.second = count_2; + data->params->trafo3D.first = trafo_3; + data->params->trafo3D.second = count_3; + data->params->placements.emplace_back(pv); + Transform3D tr3(start); + for(size_t k=0; k < count_3; ++k) { + Transform3D tr2 = tr3; + for(size_t j=0; j < count_2; ++j) { + Transform3D tr1 = tr2; + for(size_t i = 0; i < count_1; ++i) { + if ( !( i == 0 && j == 0 && k == 0 ) ) { + PlacedVolume ppv = + _addNode(m_element, entity, get_copy_number(m_element), detail::matrix::_transform(tr1)); + data->params->placements.emplace_back(ppv); + ppv.data()->params = data->params->addref(); + } + tr1 *= trafo_1; + } + tr2 *= trafo_2; + } + tr3 *= trafo_3; + } + return pv; } /// Set the volume's option value -void Volume::setOption(const string& opt) const { - m_element->SetOption(opt.c_str()); +const Volume& Volume::setOption(const string& opt) const { + if ( isValid() ) { + m_element->SetOption(opt.c_str()); + return *this; + } + throw runtime_error("dd4hep: Attempt to access invalid handle of type: PlacedVolume"); } /// Access the volume's option value @@ -832,8 +1059,6 @@ Material Volume::material() const { return Material(m_element->GetMedium()); } -#include "TROOT.h" - /// Set Visualization attributes to the volume const Volume& Volume::setVisAttributes(const VisAttr& attr) const { if ( attr.isValid() ) { diff --git a/DDG4/include/DDG4/Geant4GeometryInfo.h b/DDG4/include/DDG4/Geant4GeometryInfo.h index a449f1e94..9b5804630 100644 --- a/DDG4/include/DDG4/Geant4GeometryInfo.h +++ b/DDG4/include/DDG4/Geant4GeometryInfo.h @@ -82,6 +82,7 @@ namespace dd4hep { typedef std::map<const TGeoShape*, G4VSolid*> SolidMap; //typedef std::map<VisAttr, G4VisAttributes*> VisMap; //typedef std::map<Geant4PlacementPath, VolumeID> Geant4PathMap; + typedef std::map<const G4VPhysicalVolume*, PlacedVolume> G4PlacementMap; } /// Concreate class holding the relation information between geant4 objects and dd4hep objects. @@ -92,7 +93,22 @@ namespace dd4hep { */ class Geant4GeometryInfo : public TNamed, public detail::GeoHandlerTypes::GeometryInfo { public: - typedef std::vector<const G4VPhysicalVolume*> Geant4PlacementPath; + struct Placement { + VolumeID volumeID; + int flags; + }; + union PlacementFlags { + int value; + struct _flags { + unsigned parametrised:1; + unsigned replicated:1; + unsigned path_has_parametrised:1; + unsigned path_has_replicated:1; + } flags; + PlacementFlags() { this->value = 0; } + PlacementFlags(int v) { this->value = v; } + }; + typedef std::vector<const G4VPhysicalVolume*> Geant4PlacementPath; TGeoManager* manager = 0; Geant4GeometryMaps::IsotopeMap g4Isotopes; Geant4GeometryMaps::ElementMap g4Elements; @@ -102,6 +118,8 @@ namespace dd4hep { Geant4GeometryMaps::PlacementMap g4Placements; Geant4GeometryMaps::AssemblyMap g4AssemblyVolumes; Geant4GeometryMaps::VolumeImprintMap g4VolumeImprints; + Geant4GeometryMaps::G4PlacementMap g4Parameterised; + Geant4GeometryMaps::G4PlacementMap g4Replicated; struct PropertyVector { std::vector<double> bins; std::vector<double> values; @@ -118,7 +136,7 @@ namespace dd4hep { std::map<Region, G4Region*> g4Regions; std::map<VisAttr, G4VisAttributes*> g4Vis; std::map<LimitSet, G4UserLimits*> g4Limits; - std::map<Geant4PlacementPath, VolumeID> g4Paths; + std::map<Geant4PlacementPath, Placement> g4Paths; std::map<SensitiveDetector,std::set<const TGeoVolume*> > sensitives; std::map<Region, std::set<const TGeoVolume*> > regions; std::map<LimitSet, std::set<const TGeoVolume*> > limits; diff --git a/DDG4/include/DDG4/Geant4Helpers.h b/DDG4/include/DDG4/Geant4Helpers.h new file mode 100644 index 000000000..827e07e1d --- /dev/null +++ b/DDG4/include/DDG4/Geant4Helpers.h @@ -0,0 +1,88 @@ +//========================================================================== +// 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. +// +// Author : M.Frank +// +//========================================================================== + +#ifndef DDG4_GEANT4HELPERS_H +#define DDG4_GEANT4HELPERS_H + +/// Framework include files +#include <DD4hep/Objects.h> + +/// Geant4 include files +#include <G4ThreeVector.hh> +#include <G4Transform3D.hh> +#include <G4RotationMatrix.hh> +#include <geomdefs.hh> + +/// C/C++ include files +#include <functional> + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim { + + /// Extract G4 rotation matrix from ROOT matrices + G4RotationMatrix g4Rotation(const TGeoMatrix* matrix); + G4RotationMatrix g4Rotation(const TGeoMatrix& matrix); + G4RotationMatrix g4Rotation(const TGeoRotation* rot); + G4RotationMatrix g4Rotation(const TGeoRotation& rot); + G4RotationMatrix g4Rotation(const Rotation3D& rot); + G4RotationMatrix g4Rotation(const RotationZYX& rot); + + /// Generate parameterised placements in 1 dimension according to transformation delta + G4Transform3D + generate_placements(const G4Transform3D& start, + const G4Transform3D& delta, + std::size_t count, + const std::function<void(const G4Transform3D& delta)>& callback); + + /// Generate parameterised placements in 2 dimensions according to transformation delta + G4Transform3D + generate_placements(const G4Transform3D& start, + const G4Transform3D& delta1, + std::size_t count1, + const G4Transform3D& delta2, + std::size_t count2, + const std::function<void(const G4Transform3D& delta)>& callback); + + /// Generate parameterised placements in 3 dimensions according to transformation delta + G4Transform3D + generate_placements(const G4Transform3D& start, + const G4Transform3D& delta1, + std::size_t count1, + const G4Transform3D& delta2, + std::size_t count2, + const G4Transform3D& delta3, + std::size_t count3, + const std::function<void(const G4Transform3D& delta)>& callback); + + /// These conversions automatically apply the conversion from CM to MM! + + G4Transform3D g4Transform(const double translation[]); + G4Transform3D g4Transform(const double translation[], const double rotation[]); + + G4Transform3D g4Transform(const Transform3D& matrix); + G4Transform3D g4Transform(const Position& pos, const Rotation3D& rot); + G4Transform3D g4Transform(const Position& pos, const RotationZYX& rot); + + G4Transform3D g4Transform(const TGeoMatrix* matrix); + G4Transform3D g4Transform(const TGeoMatrix& matrix); + G4Transform3D g4Transform(const TGeoTranslation* translation, const TGeoRotation* rotation); + G4Transform3D g4Transform(const TGeoTranslation& translation, const TGeoRotation& rotation); + + std::pair<double, EAxis> extract_axis(const Transform3D& trafo); + + } // End namespace sim +} // End namespace dd4hep +#endif // DDG4_GEANT4HELPERS_H diff --git a/DDG4/include/DDG4/Geant4PlacementParameterisation.h b/DDG4/include/DDG4/Geant4PlacementParameterisation.h new file mode 100644 index 000000000..d9c45b426 --- /dev/null +++ b/DDG4/include/DDG4/Geant4PlacementParameterisation.h @@ -0,0 +1,98 @@ +//========================================================================== +// 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. +// +// Author : M.Frank +// +//========================================================================== +#ifndef DDG4_GEANT4PLACEMENTPARAMETERISATION_H +#define DDG4_GEANT4PLACEMENTPARAMETERISATION_H + +/// Geant4 include files +#include <G4PVParameterised.hh> +#include <G4VPVParameterisation.hh> + +/// C/C++ include files +#include <vector> + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim { + + /// Class to dump the records of the intrinsic Geant4 event model. + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4PlacementParameterisation : public G4VPVParameterisation { + public: + + /// Helper structure to cache dimension variabled from setup parameters + struct Dimension { + G4Transform3D delta { }; + G4ThreeVector translation { }; + std::size_t count { 0 }; + + public: + /// Initializing constructor + Dimension(const G4Transform3D& d, std::size_t c) + : delta(d), translation(d.getTranslation()), count(c) {} + /// Default constructor + Dimension() = default; + /// Move Constructor + Dimension(Dimension&&) = default; + /// Copy Constructor + Dimension(const Dimension&) = default; + /// Assignment operator + Dimension& operator=(Dimension&&) = default; + /// Assignment operator + Dimension& operator=(const Dimension&) = default; + }; + + using Parameters = PlacedVolume::Object::Parameterisation; + using Rotations = std::vector<G4RotationMatrix>; + using Translations = std::vector<G4ThreeVector>; + using Dimensions = std::vector<Dimension>; + + /// Setup parameters + Dimension m_start { }; + Dimensions m_dimensions { }; + EAxis m_axis { kUndefined }; + + /// Optimization flag for simple parameteristions + bool m_have_rotation { false }; + + /// Reference to the DD4hep placement handle + PlacedVolume m_placement; + const Parameters& m_params; + + /// Cached rotations and translations for complex parameterisations + mutable Rotations m_rotations; + mutable Translations m_translations; + + /// Callback to store resulting rotation + void operator()(const G4Transform3D& transform); + + public: + /// Initializing constructor + Geant4PlacementParameterisation(PlacedVolume pv); + /// Standard destructor + virtual ~Geant4PlacementParameterisation() = default; + /// Access Axis direction + EAxis axis() const { return m_axis; } + /// Access number of replicas + std::size_t count() const; + /// G4VPVParameterisation overload: Compute copy transformation + virtual void ComputeTransformation(const G4int copyNo, G4VPhysicalVolume* vol) const override; + }; + } // End namespace sim +} // End namespace dd4hep +#endif // DDG4_GEANT4PLACEMENTPARAMETERISATION_H diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h index 34e255284..c70e714df 100644 --- a/DDG4/include/DDG4/Geant4SensDetAction.h +++ b/DDG4/include/DDG4/Geant4SensDetAction.h @@ -14,9 +14,9 @@ #define DDG4_GEANT4SENSDETACTION_H // Framework include files -#include "DD4hep/Detector.h" -#include "DDG4/Geant4Action.h" -#include "DDG4/Geant4HitCollection.h" +#include <DD4hep/Detector.h> +#include <DDG4/Geant4Action.h> +#include <DDG4/Geant4HitCollection.h> // Geant4 include files #include <G4ThreeVector.hh> @@ -68,7 +68,7 @@ namespace dd4hep { virtual ~Geant4ActionSD(); public: /// Initialize the usage of a hit collection. Returns the collection identifier - virtual size_t defineCollection(const std::string& name) = 0; + virtual std::size_t defineCollection(const std::string& name) = 0; /// Access to the readout geometry of the sensitive detector virtual G4VReadOutGeometry* readoutGeometry() const = 0; /// This is a utility method which returns the hits collection ID @@ -239,16 +239,16 @@ namespace dd4hep { bool accept(const Geant4FastSimSpot* step) const; /// Initialize the usage of a single hit collection. Returns the collection ID - template <typename TYPE> size_t defineCollection(const std::string& coll_name); + template <typename TYPE> std::size_t defineCollection(const std::string& coll_name); /// Access HitCollection container names - const std::string& hitCollectionName(size_t which) const; + const std::string& hitCollectionName(std::size_t which) const; /// Retrieve the hits collection associated with this detector by its serial number - Geant4HitCollection* collection(size_t which); + Geant4HitCollection* collection(std::size_t which); /// Retrieve the hits collection associated with this detector by its collection identifier - Geant4HitCollection* collectionByID(size_t id); + Geant4HitCollection* collectionByID(std::size_t id); /// Define collections created by this sensitivie action object virtual void defineCollections(); @@ -370,24 +370,24 @@ namespace dd4hep { virtual void updateContext(Geant4Context* ctxt); /// Called at construction time of the sensitive detector to declare all hit collections - size_t defineCollections(Geant4ActionSD* sens_det); + std::size_t defineCollections(Geant4ActionSD* sens_det); /// Initialize the usage of a hit collection. Returns the collection identifier - size_t defineCollection(Geant4Sensitive* owner, const std::string& name, create_t func); + std::size_t defineCollection(Geant4Sensitive* owner, const std::string& name, create_t func); /// Define a named collection containing hist of a specified type - template <typename TYPE> size_t defineCollection(Geant4Sensitive* owner, const std::string& coll_name) { + template <typename TYPE> std::size_t defineCollection(Geant4Sensitive* owner, const std::string& coll_name) { return defineCollection(owner, coll_name, Geant4SensDetActionSequence::_create<TYPE>); } /// Access HitCollection container names - const std::string& hitCollectionName(size_t which) const; + const std::string& hitCollectionName(std::size_t which) const; /// Retrieve the hits collection associated with this detector by its serial number - Geant4HitCollection* collection(size_t which) const; + Geant4HitCollection* collection(std::size_t which) const; /// Retrieve the hits collection associated with this detector by its collection identifier - Geant4HitCollection* collectionByID(size_t id) const; + Geant4HitCollection* collectionByID(std::size_t id) const; /// Register begin-of-event callback template <typename T> void callAtBegin(T* p, void (T::*f)(G4HCofThisEvent*)) { @@ -493,7 +493,7 @@ namespace dd4hep { }; /// Initialize the usage of a single hit collection. Returns the collection ID - template <typename TYPE> inline size_t Geant4Sensitive::defineCollection(const std::string& coll_name) { + template <typename TYPE> inline std::size_t Geant4Sensitive::defineCollection(const std::string& coll_name) { return sequence().defineCollection<TYPE>(this, coll_name); } @@ -521,7 +521,7 @@ namespace dd4hep { std::string m_readoutName { }; /// Collection identifier - size_t m_collectionID { 0 }; + std::size_t m_collectionID { 0 }; /// User data block UserData m_userData { }; @@ -559,11 +559,11 @@ namespace dd4hep { * At the same time a VolumeID filter is injected at the front of the sensitive's * filter queue to ONLY act on volume IDs matching this criterium. */ - template <typename HIT> size_t declareReadoutFilteredCollection(); + template <typename HIT> std::size_t declareReadoutFilteredCollection(); /// Define readout specific hit collection. matching name must be present in readout structure template <typename HIT> - size_t defineReadoutCollection(const std::string collection_name); + std::size_t defineReadoutCollection(const std::string collection_name); /// Initialization overload for specialization virtual void initialize() final; @@ -586,5 +586,4 @@ namespace dd4hep { } // End namespace sim } // End namespace dd4hep - #endif // DDG4_GEANT4SENSDETACTION_H diff --git a/DDG4/include/DDG4/Geant4SensDetAction.inl b/DDG4/include/DDG4/Geant4SensDetAction.inl index 0d60d4037..235e7b14f 100644 --- a/DDG4/include/DDG4/Geant4SensDetAction.inl +++ b/DDG4/include/DDG4/Geant4SensDetAction.inl @@ -12,13 +12,13 @@ //========================================================================== /// Framework include files -#include "DDG4/Geant4SensDetAction.h" -#include "DD4hep/Detector.h" -#include "DD4hep/InstanceCount.h" -#include "DD4hep/detail/ObjectsInterna.h" +#include <DDG4/Geant4SensDetAction.h> +#include <DD4hep/Detector.h> +#include <DD4hep/InstanceCount.h> +#include <DD4hep/detail/ObjectsInterna.h> -#include "DDG4/Geant4ReadoutVolumeFilter.h" -#include "DDG4/Geant4Data.h" +#include <DDG4/Geant4ReadoutVolumeFilter.h> +#include <DDG4/Geant4Data.h> /// Namespace for the AIDA detector description toolkit namespace dd4hep { @@ -90,11 +90,11 @@ namespace dd4hep { /// Define readout specific hit collection. matching name must be present in readout structure template <typename T> template <typename HIT> - size_t Geant4SensitiveAction<T>::defineReadoutCollection(const std::string match_name) { + std::size_t Geant4SensitiveAction<T>::defineReadoutCollection(const std::string match_name) { auto ro = readout(); for(const auto& c : ro->hits ) { if ( c.name == match_name ) { - size_t coll_id = defineCollection<HIT>(match_name); + std::size_t coll_id = defineCollection<HIT>(match_name); Geant4Filter* filter = new Geant4ReadoutVolumeFilter(context(),name()+"_"+c.name,ro,c.name); adopt_front(filter); return coll_id; @@ -107,7 +107,7 @@ namespace dd4hep { /// Define readout specific hit collection with volume ID filtering template <typename T> template <typename HIT> - size_t Geant4SensitiveAction<T>::declareReadoutFilteredCollection() { + std::size_t Geant4SensitiveAction<T>::declareReadoutFilteredCollection() { if ( m_collectionName.empty() ) { return defineCollection<HIT>(readout().name()); } @@ -132,7 +132,7 @@ namespace dd4hep { } // End namespace dd4hep -#include "DD4hep/Printout.h" -#include "DDG4/Geant4TouchableHandler.h" -#include "DDG4/Geant4StepHandler.h" -#include "DDG4/Geant4VolumeManager.h" +#include <DD4hep/Printout.h> +#include <DDG4/Geant4TouchableHandler.h> +#include <DDG4/Geant4StepHandler.h> +#include <DDG4/Geant4VolumeManager.h> diff --git a/DDG4/include/DDG4/Geant4Vertex.h b/DDG4/include/DDG4/Geant4Vertex.h index b4fdd80bb..a2ee64001 100644 --- a/DDG4/include/DDG4/Geant4Vertex.h +++ b/DDG4/include/DDG4/Geant4Vertex.h @@ -14,7 +14,7 @@ #define DDG4_GEANT4VERTEX_H // Framework include files -#include "DD4hep/Memory.h" +#include <DD4hep/Memory.h> // C/C++ include files #include <set> diff --git a/DDG4/include/DDG4/Geant4VolumeManager.h b/DDG4/include/DDG4/Geant4VolumeManager.h index 7b1035d89..6cbee9517 100644 --- a/DDG4/include/DDG4/Geant4VolumeManager.h +++ b/DDG4/include/DDG4/Geant4VolumeManager.h @@ -14,9 +14,9 @@ #define DDG4_GEANT4VOLUMEMANAGER_H // Framework include files -#include "DD4hep/Detector.h" -#include "DD4hep/IDDescriptor.h" -#include "DDG4/Geant4Primitives.h" +#include <DD4hep/Detector.h> +#include <DD4hep/IDDescriptor.h> +#include <DDG4/Geant4Primitives.h> // Geant4 forward declarations class G4VTouchable; @@ -41,9 +41,6 @@ namespace dd4hep { */ class Geant4VolumeManager : public Handle<Geant4GeometryInfo> { protected: - /// Optimization flag to shortcut object checks - mutable bool m_isValid = false; - /// Check the validity of the information before accessing it. bool checkValidity() const; @@ -58,12 +55,12 @@ namespace dd4hep { Geant4VolumeManager() = default; /// Constructor to be used when reading the already parsed object Geant4VolumeManager(const Handle<Geant4GeometryInfo>& e) - : Handle<Geant4GeometryInfo>(e), m_isValid(false) { } + : Handle<Geant4GeometryInfo>(e) { } /// Constructor to be used when reading the already parsed object Geant4VolumeManager(const Geant4VolumeManager& e) = default; /// Constructor to be used when reading the already parsed object template <typename Q> Geant4VolumeManager(const Handle<Q>& e) - : Handle<Geant4GeometryInfo>(e), m_isValid(false) { + : Handle<Geant4GeometryInfo>(e) { } /// Assignment operator Geant4VolumeManager& operator=(const Geant4VolumeManager& c) = default; @@ -72,11 +69,11 @@ namespace dd4hep { std::vector<const G4VPhysicalVolume*> placementPath(const G4VTouchable* touchable, bool exception = true) const; /// Access CELLID by placement path - VolumeID volumeID(const std::vector<const G4VPhysicalVolume*>& path) const; + //VolumeID volumeID(const std::vector<const G4VPhysicalVolume*>& path) const; /// Access CELLID by Geant4 touchable object VolumeID volumeID(const G4VTouchable* touchable) const; /// Accessfully decoded volume fields by placement path - void volumeDescriptor(const std::vector<const G4VPhysicalVolume*>& path, + void volumeDescriptor(const std::vector<const G4VPhysicalVolume*>& path, std::pair<VolumeID,std::vector<std::pair<const BitFieldElement*, VolumeID> > >& volume_desc) const; /// Access fully decoded volume fields by Geant4 touchable object void volumeDescriptor(const G4VTouchable* touchable, diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp index 9633af80c..794693e6f 100644 --- a/DDG4/src/Geant4Converter.cpp +++ b/DDG4/src/Geant4Converter.cpp @@ -11,7 +11,6 @@ // //========================================================================== - // Framework include files #include <DD4hep/Detector.h> #include <DD4hep/Plugins.h> @@ -25,20 +24,21 @@ #include <DD4hep/detail/DetectorInterna.h> #include <DDG4/Geant4Field.h> +#include <DDG4/Geant4Helpers.h> #include <DDG4/Geant4Converter.h> #include <DDG4/Geant4UserLimits.h> +#include <DDG4/Geant4PlacementParameterisation.h> #include "Geant4ShapeConverter.h" // ROOT includes #include <TMath.h> #include <TROOT.h> -//#include <TColor.h> -//#include <TGeoManager.h> #include <TGeoBoolNode.h> // Geant4 include files #include <G4Version.hh> #include <G4VisAttributes.hh> +#include <G4PVParameterised.hh> #include <G4ProductionCuts.hh> #include <G4VUserRegionInformation.hh> @@ -68,9 +68,7 @@ #if G4VERSION_NUMBER >= 1040 #include <G4MaterialPropertiesIndex.hh> #endif -#if G4VERSION_NUMBER >= 1030 #include <G4ScaledSolid.hh> -#endif #include <CLHEP/Units/SystemOfUnits.h> // C/C++ include files @@ -92,18 +90,6 @@ static constexpr const char* GEANT4_TAG_IGNORE = "Geant4-ignore"; namespace { static string indent = ""; - static Double_t s_identity_rot[] = { 1., 0., 0., 0., 1., 0., 0., 0., 1. }; - struct MyTransform3D : public G4Transform3D { - MyTransform3D(double XX, double XY, double XZ, double DX, double YX, double YY, double YZ, double DY, double ZX, double ZY, - double ZZ, double DZ) - : G4Transform3D(XX, XY, XZ, DX, YX, YY, YZ, DY, ZX, ZY, ZZ, DZ) { - } - MyTransform3D(const double* t, const double* r = s_identity_rot) - : G4Transform3D(r[0],r[1],r[2],t[0]*CM_2_MM,r[3],r[4],r[5],t[1]*CM_2_MM,r[6],r[7],r[8],t[2]*CM_2_MM) { - } - MyTransform3D(Transform3D&& copy) : Transform3D(copy) {} - }; - bool is_left_handed(const TGeoMatrix* m) { const Double_t* r = m->GetRotationMatrix(); if ( r ) { @@ -131,7 +117,6 @@ namespace { } }; - pair<double,double> g4PropertyConversion(int index) { #if G4VERSION_NUMBER >= 1040 switch(index) { @@ -420,7 +405,7 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const // We need to convert the property from TGeo units to Geant4 units auto conv = g4PropertyConversion(idx); vector<double> bins(v->bins), vals(v->values); - for(size_t i=0, count=bins.size(); i<count; ++i) + for(std::size_t i=0, count=bins.size(); i<count; ++i) bins[i] *= conv.first, vals[i] *= conv.second; G4MaterialPropertyVector* vec = @@ -428,7 +413,7 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const tab->AddProperty(named->GetName(), vec); 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) + for(std::size_t i=0, count=v->bins.size(); i<count; ++i) printout(lvl, name, " Geant4: %s %8.3g [MeV] TGeo: %8.3g [GeV] Conversion: %8.3g", named->GetName(), bins[i], v->bins[i], conv.first); } @@ -556,7 +541,6 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) c solid = convertShape<TGeoTessellated>(shape); #endif else if (isa == TGeoScaledShape::Class()) { -#if G4VERSION_NUMBER >= 1030 TGeoScaledShape* sh = (TGeoScaledShape*) shape; TGeoShape* sol = sh->GetShape(); const double* vals = sh->GetScale()->GetScale(); @@ -566,9 +550,6 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) c solid = new G4ScaledSolid(sh->GetName(), g4solid, scal); else solid = new G4ReflectedSolid(g4solid->GetName()+"_refl", g4solid, scal); -#else - except("Geant4Converter","++ TGeoScaledShape are only supported by Geant4 for versions >= 10.3"); -#endif } else if ( isa == TGeoCompositeShape::Class() ) { const TGeoCompositeShape* sh = (const TGeoCompositeShape*) shape; @@ -585,10 +566,6 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) c except("Geant4Converter","++ No right Geant4 Solid present for composite shape: %s",name.c_str()); } - //specific case! - //Ellipsoid tag preparing - //if left == TGeoScaledShape AND right == TGeoBBox - // AND if TGeoScaledShape->GetShape == TGeoSphere TGeoShape* ls = boolean->GetLeftShape(); TGeoShape* rs = boolean->GetRightShape(); if (strcmp(ls->ClassName(), "TGeoScaledShape") == 0 && @@ -617,7 +594,7 @@ void* Geant4Converter::handleSolid(const string& name, const TGeoShape* shape) c } if ( matrix->IsRotation() ) { - MyTransform3D transform(matrix->GetTranslation(),matrix->GetRotationMatrix()); + G4Transform3D transform(g4Transform(matrix)); if (oper == TGeoBoolNode::kGeoSubtraction) solid = new G4SubtractionSolid(name, left, right, transform); else if (oper == TGeoBoolNode::kGeoUnion) @@ -655,7 +632,7 @@ void* Geant4Converter::handleVolume(const string& name, const TGeoVolume* volume Volume _v(volume); if ( _v.testFlagBit(Volume::VETO_SIMU) ) { printout(lvl, "Geant4Converter", "++ Volume %s not converted [Veto'ed for simulation]",volume->GetName()); - return 0; + return nullptr; } else if (volIt == info.g4Volumes.end() ) { string n = volume->GetName(); @@ -733,12 +710,12 @@ void* Geant4Converter::handleAssembly(const string& name, const TGeoNode* node) TGeoVolume* mot_vol = node->GetVolume(); PrintLevel lvl = debugVolumes ? ALWAYS : outputLevel; if ( mot_vol->IsA() != TGeoVolumeAssembly::Class() ) { - return 0; + return nullptr; } Volume _v(mot_vol); if ( _v.testFlagBit(Volume::VETO_SIMU) ) { printout(lvl, "Geant4Converter", "++ AssemblyNode %s not converted [Veto'ed for simulation]",node->GetName()); - return 0; + return nullptr; } Geant4GeometryInfo& info = data(); Geant4AssemblyVolume* g4 = info.g4AssemblyVolumes[node]; @@ -751,7 +728,7 @@ void* Geant4Converter::handleAssembly(const string& name, const TGeoNode* node) TGeoNode* dau = mot_vol->GetNode(i); TGeoVolume* dau_vol = dau->GetVolume(); TGeoMatrix* tr = dau->GetMatrix(); - MyTransform3D transform(tr->GetTranslation(),tr->IsRotation() ? tr->GetRotationMatrix() : s_identity_rot); + G4Transform3D transform(g4Transform(tr)); if ( is_left_handed(tr) ) { G4Scale3D scale; @@ -772,7 +749,7 @@ void* Geant4Converter::handleAssembly(const string& name, const TGeoNode* node) printout(FATAL, "Geant4Converter", "+++ Invalid child assembly at %s : %d parent: %s child:%s", __FILE__, __LINE__, name.c_str(), dau->GetName()); delete g4; - return 0; + return nullptr; } g4->placeAssembly(dau, (*ia).second, transform); printout(lvl, "Geant4Converter", "+++ Assembly: AddPlacedAssembly : dau:%s " @@ -811,7 +788,7 @@ void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node) if ( _v.testFlagBit(Volume::VETO_SIMU) ) { printout(lvl, "Geant4Converter", "++ Placement %s not converted [Veto'ed for simulation]",node->GetName()); - return 0; + return nullptr; } //g4 = nullptr; if ( !g4 ) { @@ -831,7 +808,7 @@ void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node) bool node_is_reflected = is_left_handed(tr); bool node_is_assembly = vol->IsA() == TGeoVolumeAssembly::Class(); bool mother_is_assembly = mot_vol ? mot_vol->IsA() == TGeoVolumeAssembly::Class() : false; - MyTransform3D transform(tr->GetTranslation(),tr->IsRotation() ? tr->GetRotationMatrix() : s_identity_rot); + G4Transform3D transform(g4Transform(tr)); Geant4GeometryMaps::VolumeMap::const_iterator volIt = info.g4Volumes.find(mot_vol); if ( mother_is_assembly ) { @@ -871,16 +848,65 @@ void* Geant4Converter::handlePlacement(const string& name, const TGeoNode* node) else if ( node != info.manager->GetTopNode() && volIt == info.g4Volumes.end() ) { throw logic_error("Geant4Converter: Invalid mother volume found!"); } + PlacedVolume pv(node); + const auto* pv_data = pv.data(); G4LogicalVolume* g4vol = info.g4Volumes[vol]; G4LogicalVolume* g4mot = info.g4Volumes[mot_vol]; - G4PhysicalVolumesPair pvPlaced = - G4ReflectionFactory::Instance()->Place(transform, // no rotation - name, // its name - g4vol, // its logical volume - g4mot, // its mother (logical) volume - false, // no boolean operations - copy, // its copy number - checkOverlaps); + G4PhysicalVolumesPair pvPlaced { nullptr, nullptr }; + + if ( pv_data && pv_data->params && (pv_data->params->flags&Volume::REPLICATED) ) { + EAxis axis = kUndefined; + double width = 0e0, offset = 0e0; + auto flags = pv_data->params->flags; + auto count = pv_data->params->trafo1D.second; + const auto& start = pv_data->params->start.Translation().Vect(); + const auto& delta = pv_data->params->trafo1D.first.Translation().Vect(); + + if ( flags&Volume::X_axis ) + { axis = kXAxis; width = delta.X(); offset = start.X(); } + else if ( flags&Volume::Y_axis ) + { axis = kYAxis; width = delta.Y(); offset = start.Y(); } + else if ( flags&Volume::Z_axis ) + { axis = kZAxis; width = delta.Z(); offset = start.Z(); } + else + except("Geant4Converter", + "++ Replication around unknown axis is not implemented. flags: %16X", flags); + pvPlaced = + G4ReflectionFactory::Instance()->Replicate(name, // its name + g4vol, // its logical volume + g4mot, // its mother (logical) volume + axis, // its replication axis + count, // Number of replicas + width, // Distanve between 2 replicas + offset); // Placement offset in axis direction + /// Update replica list to avoid additional conversions... + auto* g4pv = pvPlaced.second ? pvPlaced.second : pvPlaced.first; + for( auto& handle : pv_data->params->placements ) + info.g4Placements[handle.ptr()] = g4pv; + } + else if ( pv_data && pv_data->params ) { + auto* g4par = new Geant4PlacementParameterisation(pv); + auto* g4pv = new G4PVParameterised(name, // its name + g4vol, // its logical volume + g4mot, // its mother (logical) volume + g4par->axis(), // its replication axis + g4par->count(), // Number of replicas + g4par); // G4 parametrization + pvPlaced = { g4pv, nullptr }; + /// Update replica list to avoid additional conversions... + for( auto& handle : pv_data->params->placements ) + info.g4Placements[handle.ptr()] = g4pv; + } + else { + pvPlaced = + G4ReflectionFactory::Instance()->Place(transform, // no rotation + name, // its name + g4vol, // its logical volume + g4mot, // its mother (logical) volume + false, // no boolean operations + 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", @@ -1101,12 +1127,12 @@ void* Geant4Converter::handleMaterialProperties(TObject* mtx) const { if ( !g4 ) { PrintLevel lvl = debugMaterials ? ALWAYS : outputLevel; g4 = new Geant4GeometryInfo::PropertyVector(); - size_t rows = matrix->GetRows(); + std::size_t rows = matrix->GetRows(); g4->name = matrix->GetName(); g4->title = matrix->GetTitle(); g4->bins.reserve(rows); g4->values.reserve(rows); - for( size_t i=0; i<rows; ++i ) { + for( std::size_t i=0; i<rows; ++i ) { g4->bins.emplace_back(matrix->Get(i,0) /* *CLHEP::eV/units::eV */); g4->values.emplace_back(matrix->Get(i,1)); } @@ -1267,7 +1293,7 @@ void* Geant4Converter::handleOpticalSurface(TObject* surface) const { // We need to convert the property from TGeo units to Geant4 units auto conv = g4PropertyConversion(idx); vector<double> bins(v->bins), vals(v->values); - for(size_t i=0, count=v->bins.size(); i<count; ++i) + for(std::size_t i=0, count=v->bins.size(); i<count; ++i) bins[i] *= conv.first, vals[i] *= conv.second; G4MaterialPropertyVector* vec = new G4MaterialPropertyVector(&bins[0], &vals[0], bins.size()); tab->AddProperty(named->GetName(), vec); @@ -1275,7 +1301,7 @@ void* Geant4Converter::handleOpticalSurface(TObject* surface) const { printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter", "++ 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) + for(std::size_t i=0, count=v->bins.size(); i<count; ++i) printout(debugSurfaces ? ALWAYS : DEBUG, named->GetName(), " Geant4: %8.3g [MeV] TGeo: %8.3g [GeV] Conversion: %8.3g", bins[i], v->bins[i], conv.first); diff --git a/DDG4/src/Geant4Helpers.cpp b/DDG4/src/Geant4Helpers.cpp new file mode 100644 index 000000000..7b905b867 --- /dev/null +++ b/DDG4/src/Geant4Helpers.cpp @@ -0,0 +1,199 @@ +//========================================================================== +// 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. +// +// Author : M.Frank +// +//========================================================================== + +// Framework include files +#include <DD4hep/Printout.h> +#include <DD4hep/DD4hepUnits.h> +#include <DDG4/Geant4Helpers.h> + +#include <CLHEP/Units/SystemOfUnits.h> + +// ROOT include files +#include <TGeoMatrix.h> + +namespace { + static constexpr const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter); + + /// Overload to access protected constructor + struct MyTransform3D : public G4Transform3D { + MyTransform3D(double XX, double XY, double XZ, double DX, double YX, double YY, double YZ, double DY, double ZX, double ZY, + double ZZ, double DZ) + : G4Transform3D(XX, XY, XZ, DX, YX, YY, YZ, DY, ZX, ZY, ZZ, DZ) { + } + MyTransform3D(const double* t, const double* r) + : G4Transform3D(r[0],r[1],r[2],t[0]*CM_2_MM,r[3],r[4],r[5],t[1]*CM_2_MM,r[6],r[7],r[8],t[2]*CM_2_MM) { + } + MyTransform3D(const double* t) + : G4Transform3D(1.0, 0.0, 0.0, t[0]*CM_2_MM, 0.0, 1.0, 0.0, t[1]*CM_2_MM, 0.0, 0.0, 1.0, t[2]*CM_2_MM) { + } + MyTransform3D(Transform3D&& copy) : Transform3D(copy) {} + }; + /// Overload to access protected constructor + class MyG4RotationMatrix : public G4RotationMatrix { + public: + MyG4RotationMatrix() : G4RotationMatrix() {} + MyG4RotationMatrix(const double* r) + : G4RotationMatrix(r[0],r[1],r[2],r[3],r[4],r[5],r[6],r[7],r[8]) { + } + }; + const double identity[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0}; + +} + +G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoMatrix* matrix) { + return matrix->IsRotation() ? MyG4RotationMatrix(matrix->GetRotationMatrix()) : MyG4RotationMatrix(); +} + +G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoMatrix& matrix) { + return g4Rotation(&matrix); +} + +G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoRotation* matrix) { + return matrix->IsRotation() ? MyG4RotationMatrix(matrix->GetRotationMatrix()) : MyG4RotationMatrix(); +} + +G4RotationMatrix dd4hep::sim::g4Rotation(const TGeoRotation& matrix) { + return g4Rotation(&matrix); +} + +G4RotationMatrix dd4hep::sim::g4Rotation(const Rotation3D& rot) { + double r[9]; + rot.GetComponents(r, r+9); + return MyG4RotationMatrix(r); +} + +G4RotationMatrix dd4hep::sim::g4Rotation(const RotationZYX& rot) { + return g4Rotation(Rotation3D(rot)); +} + + +G4Transform3D dd4hep::sim::g4Transform(const double translation[]) { + return MyTransform3D(translation); +} + +G4Transform3D dd4hep::sim::g4Transform(const double translation[], const double rotation[]) { + return MyTransform3D(translation, rotation); +} + +G4Transform3D dd4hep::sim::g4Transform(const TGeoMatrix& matrix) { + return g4Transform(&matrix); +} + +G4Transform3D dd4hep::sim::g4Transform(const TGeoMatrix* matrix) { + return matrix->IsRotation() + ? g4Transform(matrix->GetTranslation(), matrix->GetRotationMatrix()) + : g4Transform(matrix->GetTranslation()); +} + +G4Transform3D dd4hep::sim::g4Transform(const TGeoTranslation& translation, const TGeoRotation& rotation) { + return g4Transform(&translation, &rotation); +} + +G4Transform3D dd4hep::sim::g4Transform(const TGeoTranslation* translation, const TGeoRotation* rotation) { + return rotation->IsRotation() + ? g4Transform(translation->GetTranslation(), rotation->GetRotationMatrix()) + : g4Transform(translation->GetTranslation()); +} + +G4Transform3D dd4hep::sim::g4Transform(const Position& pos, const Rotation3D& rot) { + double r[9], t[3] = {pos.X(), pos.Y(), pos.Z()}; + rot.GetComponents(r, r+9); + return g4Transform(t, r); +} + +G4Transform3D dd4hep::sim::g4Transform(const Position& pos, const RotationZYX& rot) { + return g4Transform(pos, Rotation3D(rot)); +} + +G4Transform3D dd4hep::sim::g4Transform(const Transform3D& matrix) { + Position pos; + Rotation3D rot; + matrix.GetDecomposition(rot, pos); + return g4Transform(pos, rot); +} + +/// Generate parameterised placements in 2 dimension according to transformation delta +G4Transform3D +dd4hep::sim::generate_placements(const G4Transform3D& start, + const G4Transform3D& delta, + std::size_t count, + const std::function<void(const G4Transform3D& delta)>& callback) +{ + G4Transform3D transform(start); + for( std::size_t i = 0; i < count; ++i ) { + callback(transform); + transform = transform * delta; + } + return transform; +} + +/// Generate parameterised placements in 2 dimensions according to transformation delta +G4Transform3D +dd4hep::sim::generate_placements(const G4Transform3D& start, + const G4Transform3D& delta1, + std::size_t count1, + const G4Transform3D& delta2, + std::size_t count2, + const std::function<void(const G4Transform3D& delta)>& callback) +{ + G4Transform3D transform2 = start; + for(std::size_t j = 0; j < count2; ++j) { + G4Transform3D transform1 = transform2; + generate_placements(transform1, delta1, count1, callback); + transform2 = transform2 * delta2; + } + return transform2; +} + +/// Generate parameterised placements in 3 dimensions according to transformation delta +G4Transform3D +dd4hep::sim::generate_placements(const G4Transform3D& start, + const G4Transform3D& delta1, + std::size_t count1, + const G4Transform3D& delta2, + std::size_t count2, + const G4Transform3D& delta3, + std::size_t count3, + const std::function<void(const G4Transform3D& delta)>& callback) +{ + G4Transform3D transform3 = start; + for(std::size_t k = 0; k < count3; ++k) { + G4Transform3D transform2 = transform3; + generate_placements(transform2, delta1, count1, delta2, count2, callback); + transform3 = transform3 * delta3; + } + return transform3; +} + +std::pair<double, EAxis> dd4hep::sim::extract_axis(const Transform3D& trafo) { + constexpr double eps = 2e-14; + G4Transform3D tr = g4Transform(trafo); + if ( fabs(tr.xx()) > (1e0-eps) && + fabs(tr.yy()) > (1e0-eps) && + fabs(tr.zz()) > (1e0-eps) ) { + if ( fabs(tr.dx()) > eps ) + return { tr.dx(), kYAxis }; + else if ( fabs(tr.dy()) > eps ) + return { tr.dy(), kYAxis }; + else if ( fabs(tr.dz()) > eps ) + return { tr.dz(), kZAxis }; + } + else if ( tr.getTranslation().mag() > eps && fabs(tr.dz()) < eps ) { + return { tr.getTranslation().rho(), kRho }; + } + else if ( fabs(tr.dz()) < eps ) { + return { tr.getRotation().phi(), kPhi }; // Is this correct ? + } + except("Geant4Converter","Invalid volume parametrization matrix. Unknown Axis!"); + return { 0e0, kUndefined }; +} diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp index 392218856..02527d492 100644 --- a/DDG4/src/Geant4InputHandling.cpp +++ b/DDG4/src/Geant4InputHandling.cpp @@ -12,19 +12,19 @@ //========================================================================== // Framework include files -#include "DDG4/Geant4InputHandling.h" -#include "DDG4/Geant4Primary.h" -#include "DDG4/Geant4Context.h" -#include "DDG4/Geant4Action.h" -#include "DDG4/Geant4PrimaryHandler.h" -#include "CLHEP/Units/SystemOfUnits.h" -#include "CLHEP/Units/PhysicalConstants.h" +#include <DDG4/Geant4InputHandling.h> +#include <DDG4/Geant4Primary.h> +#include <DDG4/Geant4Context.h> +#include <DDG4/Geant4Action.h> +#include <DDG4/Geant4PrimaryHandler.h> +#include <CLHEP/Units/SystemOfUnits.h> +#include <CLHEP/Units/PhysicalConstants.h> // Geant4 include files -#include "G4ParticleDefinition.hh" -#include "G4Event.hh" -#include "G4PrimaryVertex.hh" -#include "G4PrimaryParticle.hh" +#include <G4ParticleDefinition.hh> +#include <G4Event.hh> +#include <G4PrimaryVertex.hh> +#include <G4PrimaryParticle.hh> // C/C++ include files #include <stdexcept> diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp index 1db6105c1..a3dad3c72 100644 --- a/DDG4/src/Geant4Particle.cpp +++ b/DDG4/src/Geant4Particle.cpp @@ -12,17 +12,17 @@ //========================================================================== // Framework include files -#include "DD4hep/Printout.h" -#include "DD4hep/Primitives.h" -#include "DD4hep/InstanceCount.h" -#include "DDG4/Geant4Particle.h" -#include "TDatabasePDG.h" -#include "TParticlePDG.h" -#include "G4ParticleTable.hh" -#include "G4ParticleDefinition.hh" -#include "G4VProcess.hh" -#include "G4ChargedGeantino.hh" -#include "G4Geantino.hh" +#include <DD4hep/Printout.h> +#include <DD4hep/Primitives.h> +#include <DD4hep/InstanceCount.h> +#include <DDG4/Geant4Particle.h> +#include <TDatabasePDG.h> +#include <TParticlePDG.h> +#include <G4ParticleTable.hh> +#include <G4ParticleDefinition.hh> +#include <G4VProcess.hh> +#include <G4ChargedGeantino.hh> +#include <G4Geantino.hh> #include <sstream> #include <iostream> diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp index 4e09bb387..d99744eb3 100644 --- a/DDG4/src/Geant4ParticleHandler.cpp +++ b/DDG4/src/Geant4ParticleHandler.cpp @@ -12,27 +12,27 @@ //========================================================================== // Framework include files -#include "DD4hep/Primitives.h" -#include "DD4hep/InstanceCount.h" -#include "DDG4/Geant4StepHandler.h" -#include "DDG4/Geant4TrackHandler.h" -#include "DDG4/Geant4EventAction.h" -#include "DDG4/Geant4SensDetAction.h" -#include "DDG4/Geant4TrackingAction.h" -#include "DDG4/Geant4SteppingAction.h" -#include "DDG4/Geant4ParticleHandler.h" -#include "DDG4/Geant4UserParticleHandler.h" +#include <DD4hep/Primitives.h> +#include <DD4hep/InstanceCount.h> +#include <DDG4/Geant4StepHandler.h> +#include <DDG4/Geant4TrackHandler.h> +#include <DDG4/Geant4EventAction.h> +#include <DDG4/Geant4SensDetAction.h> +#include <DDG4/Geant4TrackingAction.h> +#include <DDG4/Geant4SteppingAction.h> +#include <DDG4/Geant4ParticleHandler.h> +#include <DDG4/Geant4UserParticleHandler.h> // Geant4 include files -#include "G4Step.hh" -#include "G4Track.hh" -#include "G4Event.hh" -#include "G4TrackStatus.hh" -#include "G4PrimaryVertex.hh" -#include "G4PrimaryParticle.hh" -#include "G4TrackingManager.hh" -#include "G4ParticleDefinition.hh" -#include "CLHEP/Units/SystemOfUnits.h" +#include <G4Step.hh> +#include <G4Track.hh> +#include <G4Event.hh> +#include <G4TrackStatus.hh> +#include <G4PrimaryVertex.hh> +#include <G4PrimaryParticle.hh> +#include <G4TrackingManager.hh> +#include <G4ParticleDefinition.hh> +#include <CLHEP/Units/SystemOfUnits.h> // C/C++ include files #include <set> diff --git a/DDG4/src/Geant4PlacementParameterisation.cpp b/DDG4/src/Geant4PlacementParameterisation.cpp new file mode 100644 index 000000000..3261c0485 --- /dev/null +++ b/DDG4/src/Geant4PlacementParameterisation.cpp @@ -0,0 +1,102 @@ +//========================================================================== +// 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. +// +// Author : M.Frank +// +//========================================================================== + +// Framework include files +#include <DD4hep/Volumes.h> +#include <DD4hep/Printout.h> +#include <DDG4/Geant4Helpers.h> +#include <DDG4/Geant4PlacementParameterisation.h> + +// Geant4 include files +#include <G4Transform3D.hh> + + /// Initializing constructor +dd4hep::sim::Geant4PlacementParameterisation::Geant4PlacementParameterisation(PlacedVolume pv) + : G4VPVParameterisation(), m_placement(pv), m_params(*pv.data()->params) +{ + auto& dim = m_dimensions; + G4Transform3D start = g4Transform(m_params.start); + G4Transform3D tr(g4Transform(m_params.trafo1D.first)); + m_start = {tr, 0UL}; + m_have_rotation = false; + dim.emplace_back(Dimension(tr, m_params.trafo1D.second)); + m_have_rotation |= dim.back().delta.getRotation().isIdentity(); + if ( m_params.trafo2D.second > 0 ) { + tr = g4Transform(m_params.trafo2D.first); + dim.emplace_back(Dimension(tr, m_params.trafo2D.second)); + m_have_rotation |= dim.back().delta.getRotation().isIdentity(); + } + if ( m_params.trafo3D.second > 0 ) { + tr = g4Transform(m_params.trafo3D.first); + dim.emplace_back(Dimension(tr, m_params.trafo3D.second)); + m_have_rotation |= dim.back().delta.getRotation().isIdentity(); + } + if ( m_have_rotation ) { + auto callback = std::bind(&Geant4PlacementParameterisation::operator(), + *this, std::placeholders::_1); + if ( dim.size() == 1 ) + generate_placements(start, + dim[0].delta, dim[0].count, callback); + else if ( dim.size() == 2 ) + generate_placements(start, + dim[0].delta, dim[0].count, + dim[1].delta, dim[1].count, callback); + else if ( dim.size() == 3 ) + generate_placements(start, + dim[0].delta, dim[0].count, + dim[1].delta, dim[1].count, + dim[2].delta, dim[2].count, callback); + } +} + +/// Access number of replicas +std::size_t dd4hep::sim::Geant4PlacementParameterisation::count() const { + std::size_t ncell = 1; + for(const auto& d : m_dimensions) ncell *= d.count; + return ncell; +} + +/// Callback to store resulting rotation +void dd4hep::sim::Geant4PlacementParameterisation::operator()(const G4Transform3D& transform) { + m_translations.emplace_back(transform.getTranslation()); + if ( this->m_have_rotation ) { + G4RotationMatrix rot = transform.inverse().getRotation(); + m_rotations.emplace_back(rot); + } +} + +/// G4VPVParameterisation overload: Callback to place sub-volumes +void dd4hep::sim::Geant4PlacementParameterisation::ComputeTransformation(const G4int copy, G4VPhysicalVolume *pv) const { + const auto& dim = m_dimensions; + std::size_t nd = dim.size(); + if ( !m_have_rotation ) { + G4ThreeVector tra = m_start.translation; + if ( nd >= 1 ) { + std::size_t dim1 = nd == 1 ? copy : (nd == 2 ? copy%dim[1].count : (nd == 3 ? copy%(dim[1].count*dim[2].count) : 0)); + tra = dim[0].translation * dim1; + } + if ( nd >= 2 ) { + std::size_t dim2 = nd == 3 ? copy%dim[2].count / dim[0].count : (nd==2 ? copy / dim[0].count : 0); + tra = tra + dim[1].translation * dim2; + } + if ( nd >= 3 ) { + std::size_t dim3 = nd == 3 ? copy / (dim[0].count*dim[1].count) : 0; + tra = tra + dim[2].translation * dim3; + } + pv->SetTranslation(tra); + return; + } + G4RotationMatrix& rot = m_rotations.at(copy); + pv->SetTranslation(m_translations.at(copy)); + pv->SetRotation(&rot); +} diff --git a/DDG4/src/Geant4SensDetAction.cpp b/DDG4/src/Geant4SensDetAction.cpp index 42d56fb12..68b1b43dd 100644 --- a/DDG4/src/Geant4SensDetAction.cpp +++ b/DDG4/src/Geant4SensDetAction.cpp @@ -12,15 +12,15 @@ //========================================================================== // Framework include files -#include "DD4hep/Printout.h" -#include "DD4hep/Primitives.h" -#include "DD4hep/InstanceCount.h" -#include "DDG4/Geant4Kernel.h" -#include "DDG4/Geant4Mapping.h" -#include "DDG4/Geant4StepHandler.h" -#include "DDG4/Geant4SensDetAction.h" -#include "DDG4/Geant4VolumeManager.h" -#include "DDG4/Geant4MonteCarloTruth.h" +#include <DD4hep/Printout.h> +#include <DD4hep/Primitives.h> +#include <DD4hep/InstanceCount.h> +#include <DDG4/Geant4Kernel.h> +#include <DDG4/Geant4Mapping.h> +#include <DDG4/Geant4StepHandler.h> +#include <DDG4/Geant4SensDetAction.h> +#include <DDG4/Geant4VolumeManager.h> +#include <DDG4/Geant4MonteCarloTruth.h> // Geant4 include files #include <G4Step.hh> @@ -190,17 +190,17 @@ Detector& Geant4Sensitive::detectorDescription() const { } /// Access HitCollection container names -const string& Geant4Sensitive::hitCollectionName(size_t which) const { +const string& Geant4Sensitive::hitCollectionName(std::size_t which) const { return sequence().hitCollectionName(which); } /// Retrieve the hits collection associated with this detector by its serial number -Geant4HitCollection* Geant4Sensitive::collection(size_t which) { +Geant4HitCollection* Geant4Sensitive::collection(std::size_t which) { return sequence().collection(which); } /// Retrieve the hits collection associated with this detector by its collection identifier -Geant4HitCollection* Geant4Sensitive::collectionByID(size_t id) { +Geant4HitCollection* Geant4Sensitive::collectionByID(std::size_t id) { return sequence().collectionByID(id); } @@ -343,14 +343,14 @@ void Geant4SensDetActionSequence::adopt(Geant4Filter* filter) { } /// Initialize the usage of a hit collection. Returns the collection identifier -size_t Geant4SensDetActionSequence::defineCollection(Geant4Sensitive* owner, const std::string& collection_name, create_t func) { +std::size_t Geant4SensDetActionSequence::defineCollection(Geant4Sensitive* owner, const std::string& collection_name, create_t func) { m_collections.emplace_back(collection_name, make_pair(owner,func)); return m_collections.size() - 1; } /// Called at construction time of the sensitive detector to declare all hit collections -size_t Geant4SensDetActionSequence::Geant4SensDetActionSequence::defineCollections(Geant4ActionSD* sens_det) { - size_t count = 0; +std::size_t Geant4SensDetActionSequence::Geant4SensDetActionSequence::defineCollections(Geant4ActionSD* sens_det) { + std::size_t count = 0; m_detector = sens_det; m_actors(&Geant4Sensitive::setDetector, sens_det); m_actors(&Geant4Sensitive::defineCollections); @@ -362,7 +362,7 @@ size_t Geant4SensDetActionSequence::Geant4SensDetActionSequence::defineCollectio } /// Access HitCollection container names -const std::string& Geant4SensDetActionSequence::hitCollectionName(size_t which) const { +const std::string& Geant4SensDetActionSequence::hitCollectionName(std::size_t which) const { if (which < m_collections.size()) { return m_collections[which].first; } @@ -372,7 +372,7 @@ const std::string& Geant4SensDetActionSequence::hitCollectionName(size_t which) } /// Retrieve the hits collection associated with this detector by its serial number -Geant4HitCollection* Geant4SensDetActionSequence::collection(size_t which) const { +Geant4HitCollection* Geant4SensDetActionSequence::collection(std::size_t which) const { if (which < m_collections.size()) { int hc_id = m_detector->GetCollectionID(which); Geant4HitCollection* c = (Geant4HitCollection*) m_hce->GetHC(hc_id); @@ -385,7 +385,7 @@ Geant4HitCollection* Geant4SensDetActionSequence::collection(size_t which) const } /// Retrieve the hits collection associated with this detector by its collection identifier -Geant4HitCollection* Geant4SensDetActionSequence::collectionByID(size_t id) const { +Geant4HitCollection* Geant4SensDetActionSequence::collectionByID(std::size_t id) const { Geant4HitCollection* c = (Geant4HitCollection*) m_hce->GetHC(id); if (c) return c; @@ -435,7 +435,7 @@ bool Geant4SensDetActionSequence::processFastSim(const Geant4FastSimSpot* spot, */ void Geant4SensDetActionSequence::begin(G4HCofThisEvent* hce) { m_hce = hce; - for (size_t count = 0; count < m_collections.size(); ++count) { + for (std::size_t count = 0; count < m_collections.size(); ++count) { const HitCollection& cr = m_collections[count]; Geant4HitCollection* col = (*cr.second.second)(name(), cr.first, cr.second.first); int id = m_detector->GetCollectionID(count); diff --git a/DDG4/src/Geant4ShapeConverter.cpp b/DDG4/src/Geant4ShapeConverter.cpp index 6f4af522d..5a5d4830f 100644 --- a/DDG4/src/Geant4ShapeConverter.cpp +++ b/DDG4/src/Geant4ShapeConverter.cpp @@ -12,43 +12,40 @@ //========================================================================== // Framework include files -#include "DD4hep/Shapes.h" -#include "DD4hep/Printout.h" -#include "DD4hep/DD4hepUnits.h" -#include "DD4hep/detail/ShapesInterna.h" +#include <DD4hep/Shapes.h> +#include <DD4hep/Printout.h> +#include <DD4hep/detail/ShapesInterna.h> #include "Geant4ShapeConverter.h" // ROOT includes -#include "TClass.h" -#include "TGeoMatrix.h" -#include "TGeoBoolNode.h" -#include "TGeoScaledShape.h" +#include <TClass.h> +#include <TGeoMatrix.h> +#include <TGeoBoolNode.h> +#include <TGeoScaledShape.h> // Geant4 include files -#include "G4Box.hh" -#include "G4Trd.hh" -#include "G4Tubs.hh" -#include "G4Trap.hh" -#include "G4Cons.hh" -#include "G4Hype.hh" -#include "G4Torus.hh" -#include "G4Sphere.hh" -#include "G4CutTubs.hh" -#include "G4Polycone.hh" -#include "G4Polyhedra.hh" -#include "G4Ellipsoid.hh" -#include "G4Paraboloid.hh" -#include "G4TwistedTubs.hh" -#include "G4GenericTrap.hh" -#include "G4ExtrudedSolid.hh" -#include "G4EllipticalTube.hh" +#include <G4Box.hh> +#include <G4Trd.hh> +#include <G4Tubs.hh> +#include <G4Trap.hh> +#include <G4Cons.hh> +#include <G4Hype.hh> +#include <G4Torus.hh> +#include <G4Sphere.hh> +#include <G4CutTubs.hh> +#include <G4Polycone.hh> +#include <G4Polyhedra.hh> +#include <G4Ellipsoid.hh> +#include <G4Paraboloid.hh> +#include <G4TwistedTubs.hh> +#include <G4GenericTrap.hh> +#include <G4ExtrudedSolid.hh> +#include <G4EllipticalTube.hh> // C/C++ include files -using namespace std; using namespace dd4hep::detail; -namespace units = dd4hep; /// Namespace for the AIDA detector description toolkit namespace dd4hep { @@ -61,9 +58,9 @@ namespace dd4hep { /// Convert a specific TGeo shape into the geant4 equivalent template <typename T> G4VSolid* convertShape(const TGeoShape* shape) { if ( shape ) { - dd4hep::except("convertShape","Unsupported shape: %s",shape->IsA()->GetName()); + except("convertShape","Unsupported shape: %s",shape->IsA()->GetName()); } - dd4hep::except("convertShape","Invalid shape conversion requested."); + except("convertShape","Invalid shape conversion requested."); return 0; } @@ -135,26 +132,26 @@ namespace dd4hep { } template <> G4VSolid* convertShape<TGeoArb8>(const TGeoShape* shape) { - vector<G4TwoVector> vertices; + std::vector<G4TwoVector> vertices; TGeoArb8* sh = (TGeoArb8*) shape; Double_t* vtx_xy = sh->GetVertices(); - for ( size_t i=0; i<8; ++i, vtx_xy +=2 ) + for ( std::size_t i=0; i<8; ++i, vtx_xy +=2 ) vertices.emplace_back(vtx_xy[0] * CM_2_MM, vtx_xy[1] * CM_2_MM); return new G4GenericTrap(sh->GetName(), sh->GetDz() * CM_2_MM, vertices); } template <> G4VSolid* convertShape<TGeoXtru>(const TGeoShape* shape) { const TGeoXtru* sh = (const TGeoXtru*) shape; - size_t nz = sh->GetNz(); - vector<G4ExtrudedSolid::ZSection> z; + std::size_t nz = sh->GetNz(); + std::vector<G4ExtrudedSolid::ZSection> z; z.reserve(nz); - for(size_t i=0; i<nz; ++i) { + for(std::size_t i=0; i<nz; ++i) { z.emplace_back(G4ExtrudedSolid::ZSection(sh->GetZ(i) * CM_2_MM, {sh->GetXOffset(i), sh->GetYOffset(i)}, sh->GetScale(i))); } - size_t np = sh->GetNvert(); - vector<G4TwoVector> polygon; + std::size_t np = sh->GetNvert(); + std::vector<G4TwoVector> polygon; polygon.reserve(np); - for(size_t i=0; i<np; ++i) { + for(std::size_t i=0; i<np; ++i) { polygon.emplace_back(sh->GetX(i) * CM_2_MM,sh->GetY(i) * CM_2_MM); } return new G4ExtrudedSolid(sh->GetName(), polygon, z); @@ -162,7 +159,7 @@ namespace dd4hep { template <> G4VSolid* convertShape<TGeoPgon>(const TGeoShape* shape) { const TGeoPgon* sh = (const TGeoPgon*) shape; - vector<double> rmin, rmax, z; + std::vector<double> rmin, rmax, z; for (Int_t i = 0; i < sh->GetNz(); ++i) { rmin.emplace_back(sh->GetRmin(i) * CM_2_MM); rmax.emplace_back(sh->GetRmax(i) * CM_2_MM); @@ -174,7 +171,7 @@ namespace dd4hep { template <> G4VSolid* convertShape<TGeoPcon>(const TGeoShape* shape) { const TGeoPcon* sh = (const TGeoPcon*) shape; - vector<double> rmin, rmax, z; + std::vector<double> rmin, rmax, z; for (Int_t i = 0; i < sh->GetNz(); ++i) { rmin.emplace_back(sh->GetRmin(i) * CM_2_MM); rmax.emplace_back(sh->GetRmax(i) * CM_2_MM); @@ -224,10 +221,10 @@ namespace dd4hep { } template <> G4VSolid* convertShape<G4GenericTrap>(const TGeoShape* shape) { - vector<G4TwoVector> vertices; + std::vector<G4TwoVector> vertices; TGeoTrap* sh = (TGeoTrap*) shape; Double_t* vtx_xy = sh->GetVertices(); - for ( size_t i=0; i<8; ++i, vtx_xy +=2 ) + for ( std::size_t i=0; i<8; ++i, vtx_xy +=2 ) vertices.emplace_back(vtx_xy[0] * CM_2_MM, vtx_xy[1] * CM_2_MM); return new G4GenericTrap(sh->GetName(), sh->GetDz() * CM_2_MM, vertices); } @@ -235,9 +232,9 @@ namespace dd4hep { } // End namespace dd4hep #if ROOT_VERSION_CODE > ROOT_VERSION(6,21,0) -#include "G4TessellatedSolid.hh" -#include "G4TriangularFacet.hh" -#include "G4QuadrangularFacet.hh" +#include <G4TessellatedSolid.hh> +#include <G4TriangularFacet.hh> +#include <G4QuadrangularFacet.hh> /// Namespace for the AIDA detector description toolkit namespace dd4hep { @@ -282,4 +279,5 @@ namespace dd4hep { } // End namespace sim } // End namespace dd4hep -#endif +#endif // ROOT_VERSION_CODE > ROOT_VERSION(6,21,0) + diff --git a/DDG4/src/Geant4VolumeManager.cpp b/DDG4/src/Geant4VolumeManager.cpp index a9670cc15..7449c5a13 100644 --- a/DDG4/src/Geant4VolumeManager.cpp +++ b/DDG4/src/Geant4VolumeManager.cpp @@ -12,18 +12,18 @@ //========================================================================== // Framework include files -#include "DD4hep/Printout.h" -#include "DD4hep/Volumes.h" -#include "DD4hep/DetElement.h" -#include "DD4hep/DetectorTools.h" -#include "DDG4/Geant4VolumeManager.h" -#include "DDG4/Geant4TouchableHandler.h" -#include "DDG4/Geant4Mapping.h" +#include <DD4hep/Printout.h> +#include <DD4hep/Volumes.h> +#include <DD4hep/DetElement.h> +#include <DD4hep/DetectorTools.h> +#include <DDG4/Geant4VolumeManager.h> +#include <DDG4/Geant4TouchableHandler.h> +#include <DDG4/Geant4Mapping.h> // Geant4 include files -#include "G4VTouchable.hh" -#include "G4LogicalVolume.hh" -#include "G4VPhysicalVolume.hh" +#include <G4VTouchable.hh> +#include <G4LogicalVolume.hh> +#include <G4VPhysicalVolume.hh> // C/C++ include files #include <sstream> @@ -35,7 +35,7 @@ using namespace dd4hep::sim; using namespace dd4hep; using namespace std; -#include "DDG4/Geant4AssemblyVolume.h" +#include <DDG4/Geant4AssemblyVolume.h> typedef pair<VolumeID,vector<pair<const BitFieldElement*, VolumeID> > > VolIDDescriptor; namespace { @@ -72,6 +72,13 @@ namespace { } printout(WARNING, "Geant4VolumeManager", "++ Detector element %s of type %s has no placement.", de.name(), de.type().c_str()); } + /// Needed to compute the cellID of parameterized volumes + for( const auto& pv : m_geo.g4Placements ) { + if ( pv.second->IsParameterised() ) + m_geo.g4Parameterised[pv.second] = pv.first; + if ( pv.second->IsReplicated() ) + m_geo.g4Replicated[pv.second] = pv.first; + } } /// Scan a single physical volume and look for sensitive elements below @@ -103,7 +110,7 @@ namespace { chain.pop_back(); } - void add_entry(SensitiveDetector sd, const TGeoNode* /* n */, const PlacedVolume::VolIDs& ids, const Chain& nodes) { + void add_entry(SensitiveDetector sd, const TGeoNode* n, const PlacedVolume::VolIDs& ids, const Chain& nodes) { Chain control; const TGeoNode* node; Volume vol; @@ -126,9 +133,17 @@ namespace { node = *(k); PlacementMap::const_iterator g4pit = m_geo.g4Placements.find(node); if (g4pit != m_geo.g4Placements.end()) { - path.emplace_back((*g4pit).second); + G4VPhysicalVolume* phys = (*g4pit).second; + if ( phys->IsParameterised() ) { + PlacedVolume pv(n); + PlacedVolumeExtension* ext = pv.data(); + if ( nullptr == ext->params->field ) { + ext->params->field = iddesc.field(ext->volIDs.at(0).first); + } + } + path.emplace_back(phys); printout(print_chain, "Geant4VolumeManager", "+++ Chain: Node OK: %s [%s]", - node->GetName(), (*g4pit).second->GetName().c_str()); + node->GetName(), phys->GetName().c_str()); continue; } control.insert(control.begin(),node); @@ -136,7 +151,6 @@ namespace { VolumeImprintMap::const_iterator iVolImp = m_geo.g4VolumeImprints.find(vol); if ( iVolImp != m_geo.g4VolumeImprints.end() ) { const Imprints& imprints = (*iVolImp).second; - //size_t len = kend-k; for(const auto& imp : imprints ) { const VolumeChain& c = imp.first; if ( c.size() <= control.size() && control == c ) { @@ -156,11 +170,22 @@ namespace { path.erase(path.begin()+path.size()-1); printout(print_res, "Geant4VolumeManager", "+++ Map %016X to Geant4 Path:%s", (void*)code, Geant4GeometryInfo::placementPath(path).c_str()); - if (m_geo.g4Paths.find(path) == m_geo.g4Paths.end()) { - m_geo.g4Paths[path] = code; + if ( m_geo.g4Paths.find(path) == m_geo.g4Paths.end() ) { + Geant4GeometryInfo::PlacementFlags opt; + for(const auto* phys : path) { + opt.flags.path_has_parametrised = phys->IsParameterised() ? 1 : 0; + opt.flags.path_has_replicated = phys->IsReplicated() ? 1 : 0; + } + opt.flags.parametrised = path.front()->IsParameterised() ? 1 : 0; + opt.flags.replicated = path.front()->IsReplicated() ? 1 : 0; + m_geo.g4Paths[path] = { code, opt.value }; m_entries.emplace(code,path); return; } + /// This is a normal case for parametrized volumes and no error + if ( !path.empty() && (path.front()->IsParameterised() || path.front()->IsReplicated()) ) { + return; + } printout(ERROR, "Geant4VolumeManager", "populate: Severe error: Duplicated Geant4 path!!!! %s %s", " [THIS SHOULD NEVER HAPPEN]",Geant4GeometryInfo::placementPath(path).c_str()); goto Err; @@ -169,6 +194,12 @@ namespace { int(control.size()),detail::tools::placementPath(control,true).c_str()); goto Err; } + else { + /// This is a normal case for parametrized volumes and no error + if ( !path.empty() && (path.front()->IsParameterised() || path.front()->IsReplicated()) ) { + return; + } + } printout(ERROR, "Geant4VolumeManager", "populate: Severe error: Duplicated Volume entry: 0x%X" " [THIS SHOULD NEVER HAPPEN]", code); @@ -187,7 +218,7 @@ namespace { /// Initializing constructor. The tree will automatically be built if possible Geant4VolumeManager::Geant4VolumeManager(const Detector& description, Geant4GeometryInfo* info) - : Handle<Geant4GeometryInfo>(info), m_isValid(false) { + : Handle<Geant4GeometryInfo>(info) { if (info && info->valid && info->g4Paths.empty()) { Populator p(description, *info); p.populate(description.world()); @@ -205,26 +236,24 @@ Geant4VolumeManager::placementPath(const G4VTouchable* touchable, bool exception /// Check the validity of the information before accessing it. bool Geant4VolumeManager::checkValidity() const { - if (m_isValid) { - return true; - } - else if (!isValid()) { + if (!isValid()) { throw runtime_error(format("Geant4VolumeManager", "Attempt to use invalid Geant4 volume manager [Invalid-Handle]")); } else if (!ptr()->valid) { throw runtime_error(format("Geant4VolumeManager", "Attempt to use invalid Geant4 geometry info [Invalid-Info]")); } - m_isValid = true; - return m_isValid; + return true; } +#if 0 /// Access CELLID by placement path VolumeID Geant4VolumeManager::volumeID(const vector<const G4VPhysicalVolume*>& path) const { if (!path.empty() && checkValidity()) { const auto& mapping = ptr()->g4Paths; auto i = mapping.find(path); - if (i != mapping.end()) - return (*i).second; + if ( i != mapping.end() ) { + return (*i).second.first; + } if (!path[0]) return InvalidPath; else if (!path[0]->GetLogicalVolume()->GetSensitiveDetector()) @@ -234,11 +263,60 @@ VolumeID Geant4VolumeManager::volumeID(const vector<const G4VPhysicalVolume*>& p Geant4GeometryInfo::placementPath(path).c_str()); return NonExisting; } +#endif /// Access CELLID by Geant4 touchable object VolumeID Geant4VolumeManager::volumeID(const G4VTouchable* touchable) const { Geant4TouchableHandler handler(touchable); - return volumeID(handler.placementPath()); + vector<const G4VPhysicalVolume*> path = handler.placementPath(); + if (!path.empty() && checkValidity()) { + const auto& mapping = ptr()->g4Paths; + auto i = mapping.find(path); + if ( i != mapping.end() ) { + const auto& e = (*i).second; + /// No parametrization or replication. + if ( e.flags == 0 ) { + return e.volumeID; + } + VolumeID volid = e.volumeID; + const auto& paramterised = ptr()->g4Parameterised; + const auto& replicated = ptr()->g4Replicated; + /// This is incredibly slow .... but what can I do ? Need a better idea. + for ( std::size_t j=0; j < path.size(); ++j ) { + const auto* phys = path[j]; + if ( phys->IsParameterised() ) { + int copy_no = touchable->GetCopyNumber(j); + const auto it = paramterised.find(phys); + if ( it != paramterised.end() ) { + //printout(INFO,"Geant4VolumeManager", + // "Copy number: %ld <--> %ld", copy_no, long(phys->GetCopyNo())); + const auto* field = (*it).second.data()->params->field; + volid |= IDDescriptor::encode(field, copy_no); + continue; + } + except("Geant4VolumeManager","Error Geant4VolumeManager::volumeID(const G4VTouchable* touchable)"); + } + else if ( phys->IsReplicated() ) { + int copy_no = touchable->GetCopyNumber(j); + const auto it = replicated.find(phys); + if ( it != replicated.end() ) { + const auto* field = (*it).second.data()->params->field; + volid |= IDDescriptor::encode(field, copy_no); + continue; + } + except("Geant4VolumeManager","Error Geant4VolumeManager::volumeID(const G4VTouchable* touchable)"); + } + } + return volid; + } + if (!path[0]) + return InvalidPath; + else if (!path[0]->GetLogicalVolume()->GetSensitiveDetector()) + return Insensitive; + } + printout(INFO, "Geant4VolumeManager","+++ Bad volume Geant4 Path: %s", + Geant4GeometryInfo::placementPath(path).c_str()); + return NonExisting; } /// Accessfully decoded volume fields by placement path @@ -247,17 +325,17 @@ void Geant4VolumeManager::volumeDescriptor(const vector<const G4VPhysicalVolume* { vol_desc.second.clear(); vol_desc.first = NonExisting; - if (!path.empty() && checkValidity()) { + if ( !path.empty() && checkValidity() ) { const auto& mapping = ptr()->g4Paths; auto i = mapping.find(path); if (i != mapping.end()) { - VolumeID vid = (*i).second; + VolumeID vid = (*i).second.volumeID; G4LogicalVolume* lvol = path[0]->GetLogicalVolume(); - if (lvol->GetSensitiveDetector()) { + if ( lvol->GetSensitiveDetector() ) { const G4VPhysicalVolume* node = path[0]; const PlacementMap& pm = ptr()->g4Placements; for (PlacementMap::const_iterator ipm = pm.begin(); ipm != pm.end(); ++ipm) { - if ((*ipm).second == node) { + if ( (*ipm).second == node ) { PlacedVolume pv = (*ipm).first; SensitiveDetector sd = pv.volume().sensitiveDetector(); IDDescriptor dsc = sd.readout().idSpec(); @@ -270,9 +348,9 @@ void Geant4VolumeManager::volumeDescriptor(const vector<const G4VPhysicalVolume* vol_desc.first = Insensitive; return; } - if (!path[0]) + if ( !path[0] ) vol_desc.first = InvalidPath; - else if (!path[0]->GetLogicalVolume()->GetSensitiveDetector()) + else if ( !path[0]->GetLogicalVolume()->GetSensitiveDetector() ) vol_desc.first = Insensitive; else vol_desc.first = NonExisting; diff --git a/DDG4/tpython/DDPython.cpp b/DDG4/tpython/DDPython.cpp index 098376f3f..7c02786ec 100644 --- a/DDG4/tpython/DDPython.cpp +++ b/DDG4/tpython/DDPython.cpp @@ -142,7 +142,9 @@ DDPython::DDPython() : context(0) { bool inited = ::Py_IsInitialized(); if ( !inited ) { ::Py_Initialize(); +#if PY_MAJOR_VERSION <=2 || (PY_MAJOR_VERSION == 3 && PY_MINOR_VERSION < 7) ::PyEval_InitThreads(); +#endif } else { _refCount += 1000; // Ensure we do not call Py_Finalize()! @@ -310,8 +312,10 @@ void DDPython::afterFork() const { #else ::PyOS_AfterFork_Child(); #endif +#if PY_MAJOR_VERSION <=2 || (PY_MAJOR_VERSION == 3 && PY_MINOR_VERSION < 7) ::PyEval_InitThreads(); ::PyEval_ReleaseLock(); +#endif } } diff --git a/examples/ClientTests/compact/ParamVolume1D.xml b/examples/ClientTests/compact/ParamVolume1D.xml new file mode 100644 index 000000000..7d8937cbe --- /dev/null +++ b/examples/ClientTests/compact/ParamVolume1D.xml @@ -0,0 +1,162 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd> + <!-- #========================================================================== + # 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. + # + #========================================================================== + --> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/materials.xml"/> + </includes> + + <define> + <constant name="world_size" value="30*m"/> + <constant name="world_x" value="world_size"/> + <constant name="world_y" value="world_size"/> + <constant name="world_z" value="world_size"/> + </define> + + <display> + <vis name="Invisible" showDaughters="false" visible="false"/> + <vis name="InvisibleWithChildren" showDaughters="true" visible="false"/> + <vis name="VisibleRed" r="1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="VisibleBlue" r="0.0" g="0.0" b="1.0" showDaughters="false" visible="true"/> + <vis name="VisibleYellow" r="1.0" g="1.0" b="0.0" showDaughters="false" visible="true"/> + <vis name="VisibleViolet" r="1.0" g="0.0" b="1.0" showDaughters="false" visible="true"/> + <vis name="VisibleGreen" alpha="0.1" r="0.0" g="1.0" b="0.0" drawingStyle="solid" lineStyle="solid" showDaughters="true" visible="true"/> + </display> + + <limits> + <limitset name="param_limits"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + </limitset> + </limits> + + <detectors> + <detector id="1" name="Param1D_1" type="DD4hep_ParamVolume" vis="VisibleGreen" limits="param_limits"> + <box x="20*cm" y="15*cm" z="150*cm" material="Air"/> + <param x="15*cm" y="10*cm" z="2*cm" material="Iron" vis="VisibleRed" limits="param_limits"> + <transformation> + <dim_x repeat="20"> + <position x="0*cm" y="0*cm" z="10*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_x> + </transformation> + <start> + <position x="0*cm" y="0*cm" z="-120*cm"/> + <rotation x="0" y="0" z="0"/> + </start> + </param> + <position x="0" y="0" z="175*cm"/> + <rotation x="0" y="0" z="0"/> + </detector> + <detector id="2" name="Param1D_2" type="DD4hep_ParamVolume" vis="VisibleGreen" limits="param_limits"> + <box x="20*cm" y="100*cm" z="15*cm" material="Air"/> + <param x="15*cm" y="2*cm" z="10*cm" material="Iron" vis="VisibleBlue" limits="param_limits"> + <transformation> + <dim_x repeat="12"> + <position x="0*cm" y="10*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_x> + </transformation> + <start> + <position x="0*cm" y="-50*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </start> + </param> + <position x="0" y="130*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </detector> + <detector id="3" name="Param1D_3" type="DD4hep_ParamVolume" vis="VisibleGreen" limits="param_limits"> + <box x="100*cm" y="20*cm" z="15*cm" material="Air"/> + <param x="2*cm" y="15*cm" z="10*cm" material="Iron" vis="VisibleYellow" limits="param_limits"> + <transformation> + <dim_x repeat="9"> + <position x="10*cm" y="0*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_x> + </transformation> + <start> + <position x="0*cm" y="0*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </start> + </param> + <position x="130*cm" y="0*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </detector> + + <detector id="4" name="Param1D_4" type="DD4hep_ParamVolume" vis="VisibleGreen" limits="param_limits"> + <box x="100*cm" y="20*cm" z="15*cm" material="Air"/> + <param x="2*cm" y="15*cm" z="10*cm" material="Iron" vis="VisibleViolet" limits="param_limits"> + <transformation> + <dim_x repeat="8"> + <position x="-10*cm" y="0*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_x> + </transformation> + <start> + <position x="0*cm" y="0*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </start> + </param> + <position x="-130*cm" y="0*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </detector> + + <detector id="5" name="Param1D_5" type="DD4hep_ParamVolume" vis="VisibleGreen" limits="param_limits"> + <box x="100*cm" y="100*cm" z="140*cm" material="Air"/> + <param x="2*cm" y="15*cm" z="10*cm" material="Iron" vis="VisibleViolet" limits="param_limits"> + <transformation> + <dim_x repeat="100"> + <position x="10*cm" y="25*cm" z="2.5*cm"/> + <rotation x="pi/8" y="0" z="0"/> + </dim_x> + </transformation> + <start> + <position x="40*cm" y="-40*cm" z="-130*cm"/> + <rotation x="0" y="0" z="0"/> + </start> + </param> + <position x="150*cm" y="150*cm" z="170*cm"/> + <rotation x="0" y="pi/4*0" z="0"/> + </detector> + + </detectors> + + <readouts> + <readout name="Hits1"> + <segmentation type="CartesianGridXY" grid_size_x="1*cm" grid_size_y="1*cm"/> + <id>system:6,x:12:-6,y:24:-6</id> + </readout> + <readout name="Hits2"> + <segmentation type="CartesianGridXY" grid_size_x="1*cm" grid_size_y="1*cm"/> + <id>system:6,x:12:-6,y:24:-6</id> + </readout> + <readout name="Hits3"> + <segmentation type="CartesianGridXY" grid_size_x="1*cm" grid_size_y="1*cm"/> + <id>system:6,x:12:-6,y:24:-6</id> + </readout> + <readout name="Hits4"> + <segmentation type="CartesianGridXY" grid_size_x="1*cm" grid_size_y="1*cm"/> + <id>system:6,side:2,x:12:-6,y:24:-6</id> + </readout> + </readouts> + + <fields> + <field name="GlobalSolenoid" type="solenoid" + inner_field="5.0*tesla" + outer_field="-1.5*tesla" + zmax="2*m" + outer_radius="3*m"> + </field> + </fields> + +</lccdd> diff --git a/examples/ClientTests/compact/ParamVolume2D.xml b/examples/ClientTests/compact/ParamVolume2D.xml new file mode 100644 index 000000000..037aab00e --- /dev/null +++ b/examples/ClientTests/compact/ParamVolume2D.xml @@ -0,0 +1,83 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd> + <!-- #========================================================================== + # 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. + # + #========================================================================== + --> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/materials.xml"/> + </includes> + + <define> + <constant name="world_size" value="30*m"/> + <constant name="world_x" value="world_size"/> + <constant name="world_y" value="world_size"/> + <constant name="world_z" value="world_size"/> + </define> + + <display> + <vis name="Invisible" showDaughters="false" visible="false"/> + <vis name="InvisibleWithChildren" showDaughters="true" visible="false"/> + <vis name="VisibleRed" r="1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="VisibleBlue" r="0.0" g="0.0" b="1.0" showDaughters="false" visible="true"/> + <vis name="VisibleYellow" r="1.0" g="1.0" b="0.0" showDaughters="false" visible="true"/> + <vis name="VisibleViolet" r="1.0" g="0.0" b="1.0" showDaughters="false" visible="true"/> + <vis name="VisibleGreen" alpha="0.1" r="0.0" g="1.0" b="0.0" drawingStyle="solid" lineStyle="solid" showDaughters="true" visible="true"/> + </display> + + <limits> + <limitset name="param_limits"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + </limitset> + </limits> + + <detectors> + <detector id="1" name="Param2D" type="DD4hep_ParamVolume" vis="VisibleGreen" limits="param_limits"> + <box x="20*cm" y="120*cm" z="120*cm" material="Air"/> + <param x="2*cm" y="2*cm" z="2*cm" material="Iron" vis="VisibleRed" limits="param_limits"> + <transformation> + <dim_x repeat="20"> + <position x="0*cm" y="0*cm" z="10*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_x> + <dim_y repeat="20"> + <position x="0*cm" y="10*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_y> + </transformation> + <start> + <position x="0*cm" y="-100*cm" z="-100*cm"/> + <rotation x="0" y="0" z="0"/> + </start> + </param> + <position x="0" y="0" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </detector> + </detectors> + + <readouts> + <readout name="Hits1"> + <segmentation type="CartesianGridXY" grid_size_x="1*cm" grid_size_y="1*cm"/> + <id>system:6,x:12:-6,y:24:-6</id> + </readout> + </readouts> + + <fields> + <field name="GlobalSolenoid" type="solenoid" + inner_field="5.0*tesla" + outer_field="-1.5*tesla" + zmax="2*m" + outer_radius="3*m"> + </field> + </fields> + +</lccdd> diff --git a/examples/ClientTests/compact/ParamVolume3D.xml b/examples/ClientTests/compact/ParamVolume3D.xml new file mode 100644 index 000000000..a9229f80c --- /dev/null +++ b/examples/ClientTests/compact/ParamVolume3D.xml @@ -0,0 +1,95 @@ +<?xml version="1.0" encoding="UTF-8"?> +<lccdd> + <!-- #========================================================================== + # 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. + # + #========================================================================== + --> + + <includes> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> + <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/materials.xml"/> + </includes> + + <define> + <constant name="world_size" value="30*m"/> + <constant name="world_x" value="world_size"/> + <constant name="world_y" value="world_size"/> + <constant name="world_z" value="world_size"/> + </define> + + <display> + <vis name="Invisible" showDaughters="false" visible="false"/> + <vis name="InvisibleWithChildren" showDaughters="true" visible="false"/> + <vis name="VisibleRed" r="1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/> + <vis name="VisibleBlue" r="0.0" g="0.0" b="1.0" showDaughters="false" visible="true"/> + <vis name="VisibleYellow" r="1.0" g="1.0" b="0.0" showDaughters="false" visible="true"/> + <vis name="VisibleViolet" r="1.0" g="0.0" b="1.0" showDaughters="false" visible="true"/> + <vis name="VisibleGreen" alpha="0.1" r="0.0" g="1.0" b="0.0" drawingStyle="solid" lineStyle="solid" showDaughters="true" visible="true"/> + </display> + + <limits> + <limitset name="param_limits"> + <limit name="step_length_max" particles="*" value="5.0" unit="mm" /> + </limitset> + </limits> + + <detectors> + <detector id="1" name="Param3D_1" type="DD4hep_ParamVolume" vis="VisibleGreen" readout="Hits3D" limits="param_limits"> + <!-- Mother volume dimensions --> + <box x="120*cm" y="120*cm" z="120*cm" material="Air"/> + <!-- Mother volume placement in the world --> + <position x="-60*cm" y="-60*cm" z="-60*cm"/> + <rotation x="0" y="0" z="0"/> + + <!-- Parameterisation arguments --> + <param x="2*cm" y="2*cm" z="2*cm" material="Iron" vis="VisibleRed" limits="param_limits"> + <start> + <!-- Parameterisation start position in the mother volume frame --> + <position x="-100*cm" y="-100*cm" z="-100*cm"/> + <rotation x="0" y="0" z="0"/> + </start> + <transformation> + <!-- Parameterisation in dimension 1 with 20 placements --> + <dim_x repeat="20"> + <position x="0*cm" y="0*cm" z="10*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_x> + <!-- Parameterisation in dimension 2 with 20 placements --> + <dim_y repeat="20"> + <position x="0*cm" y="10*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_y> + <!-- Parameterisation in dimension 3 with 20 placements --> + <dim_z repeat="20"> + <position x="10*cm" y="0*cm" z="0*cm"/> + <rotation x="0" y="0" z="0"/> + </dim_z> + </transformation> + </param> + </detector> + </detectors> + + <readouts> + <readout name="Hits3D"> + <segmentation type="CartesianGridXY" grid_size_x="1*mm" grid_size_y="1*mm"/> + <id>system:6,volume:16,x:32:-6,y:48:-6</id> + </readout> + </readouts> + + <fields> + <field name="GlobalSolenoid" type="solenoid" + inner_field="5.0*tesla" + outer_field="-1.5*tesla" + zmax="2*m" + outer_radius="3*m"> + </field> + </fields> + +</lccdd> diff --git a/examples/ClientTests/scripts/ParamVolume.mac b/examples/ClientTests/scripts/ParamVolume.mac new file mode 100644 index 000000000..04d245872 --- /dev/null +++ b/examples/ClientTests/scripts/ParamVolume.mac @@ -0,0 +1,11 @@ +# +#/vis/open OGLSX +# +/vis/open OGLIX +/vis/scene/create +/vis/viewer/set/background white +/vis/viewer/set/autoRefresh true +/vis/sceneHandler/attach + +/vis/viewer/set/viewpointVector 1 1 1 +/vis/scene/add/axes 0 0 0 3 m diff --git a/examples/ClientTests/scripts/ParamVolume.py b/examples/ClientTests/scripts/ParamVolume.py new file mode 100644 index 000000000..63ccbff56 --- /dev/null +++ b/examples/ClientTests/scripts/ParamVolume.py @@ -0,0 +1,124 @@ +# ========================================================================== +# 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 logging +# +logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO) +logger = logging.getLogger(__name__) +# +# +""" + + dd4hep simulation example setup using the python configuration + + @author M.Frank + @version 1.0 + +""" + + +def run(): + import os + import sys + import time + import DDG4 + args = DDG4.CommandLine() + install_dir = os.environ['DD4hepExamplesINSTALL'] + install_dir = install_dir + "/examples/ClientTests/compact" + if args.help: + logger.info(""" + python <dir>/ParamValue.py -option [-option] + -geometry <geometry file name> File is expected in the examples + install area: + """ + install_dir + """ + -vis Enable visualization + -macro Pass G4 macro file to UI executive + -batch Run in batch mode for unit testing + -event <number> Run geant4 for specified number of events + (batch mode only) + """) + sys.exit(0) + + from DDG4 import OutputLevel as Output + from g4units import GeV, MeV, mm + + kernel = DDG4.Kernel() + kernel.loadGeometry(str("file:" + install_dir + os.sep + args.geometry)) + + DDG4.importConstants(kernel.detectorDescription(), debug=False) + geant4 = DDG4.Geant4(kernel, calo='Geant4CalorimeterAction') + geant4.printDetectors() + + # Configure UI + if args.macro: + ui = geant4.setupCshUI(macro=args.macro, vis=args.vis) + else: + ui = geant4.setupCshUI(vis=args.vis) + if args.batch: + ui.Commands = ['/run/beamOn ' + str(args.events), '/ddg4/UI/terminate'] + + if args.sensitive: + geant4.setupDetectors() + + # 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.DebugElements = False + act.DebugVolumes = True + act.DebugShapes = True + seq, act = geant4.addDetectorConstruction("Geant4DetectorSensitivesConstruction/ConstructSD") + + # Configure I/O + geant4.setupROOTOutput('RootOutput', 'ParamVolume1D_' + time.strftime('%Y-%m-%d_%H-%M')) + + # Setup particle gun + gun = geant4.setupGun("Gun", particle='e+', energy=20 * GeV, multiplicity=1) + gun.OutputLevel = generator_output_level + + # And handle the simulation particles. + part = DDG4.GeneratorAction(kernel, "Geant4ParticleHandler/ParticleHandler") + kernel.generatorAction().adopt(part) + part.SaveProcesses = ['Decay'] + part.MinimalKineticEnergy = 100 * MeV + part.OutputLevel = Output.INFO # generator_output_level + part.enableUI() + user = DDG4.Action(kernel, "Geant4TCUserParticleHandler/UserParticleHandler") + user.TrackingVolume_Zmax = 1 * mm + user.TrackingVolume_Rmax = 1 * mm + user.enableUI() + part.adopt(user) + + # Now build the physics list: + phys = geant4.setupPhysics('QGSP_BERT') + ph = DDG4.PhysicsList(kernel, str('Geant4PhysicsList/Myphysics')) + ph.addParticleConstructor(str('G4Geantino')) + ph.addParticleConstructor(str('G4BosonConstructor')) + ph.enableUI() + phys.adopt(ph) + phys.dump() + + geant4.execute() + + +if __name__ == "__main__": + run() diff --git a/examples/ClientTests/scripts/SiliconBlockGFlash.py b/examples/ClientTests/scripts/SiliconBlockGFlash.py index 133b30551..e1d2cccd9 100644 --- a/examples/ClientTests/scripts/SiliconBlockGFlash.py +++ b/examples/ClientTests/scripts/SiliconBlockGFlash.py @@ -11,11 +11,10 @@ # # from __future__ import absolute_import, unicode_literals -import os -import time -import DDG4 -from DDG4 import OutputLevel as Output -from g4units import GeV, MeV, m +import logging +# +logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO) +logger = logging.getLogger(__name__) # # """ @@ -34,17 +33,38 @@ from g4units import GeV, MeV, m def run(): + import os + import time + import DDG4 args = DDG4.CommandLine() + if args.help: + import sys + logger.info(""" + python <dir>/ParamValue.py -option [-option] + -geometry <geometry file> Geometry with full path + -vis Enable visualization + -macro Pass G4 macro file to UI executive + -batch Run in batch mode for unit testing + -event <number> Run geant4 for specified number of events + (batch mode only) + """) + sys.exit(0) + + from DDG4 import OutputLevel as Output + from g4units import GeV, MeV, m kernel = DDG4.Kernel() - install_dir = os.environ['DD4hepExamplesINSTALL'] - kernel.loadGeometry(str("file:" + install_dir + "/examples/ClientTests/compact/SiliconBlock.xml")) + if args.geometry: + kernel.loadGeometry(str("file:" + args.geometry)) + else: + install_dir = os.environ['DD4hepExamplesINSTALL'] + kernel.loadGeometry(str("file:" + install_dir + "/examples/ClientTests/compact/SiliconBlock.xml")) DDG4.importConstants(kernel.detectorDescription(), debug=False) geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction', calo='Geant4CalorimeterAction') geant4.printDetectors() # Configure UI if args.macro: - ui = geant4.setupCshUI(macro=args.macro) + ui = geant4.setupCshUI(macro=args.macro, vis=args.vis) else: ui = geant4.setupCshUI() if args.batch: diff --git a/examples/ClientTests/src/ParamVolume_geo.cpp b/examples/ClientTests/src/ParamVolume_geo.cpp new file mode 100644 index 000000000..123936f38 --- /dev/null +++ b/examples/ClientTests/src/ParamVolume_geo.cpp @@ -0,0 +1,99 @@ +//========================================================================== +// 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. +// +// Author : M.Frank +// +//========================================================================== + +// Framework includes +#include "DD4hep/DetFactoryHelper.h" + +using namespace dd4hep; +using namespace dd4hep::detail; + +namespace { + Transform3D get_trafo(xml_dim_t elt) { + xml_dim_t x_pos = elt.position(); + xml_dim_t x_rot = elt.rotation(); + return Transform3D(RotationZYX(x_rot.x(), x_rot.y(), x_rot.z()), + Position(x_pos.x(), x_pos.y(), x_pos.z())); + } +} + +static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector sens) { + // XML detector object: DDCore/XML/XMLDetector.h + xml_dim_t x_det = e; + //Create the DetElement for dd4hep + DetElement det(x_det.nameStr(),x_det.id()); + + // XML dimension object: DDCore/XML/XMLDimension.h + xml_dim_t x_box(x_det.child(_U(box))); + Volume envelope_vol(x_det.nameStr()+"_envelope", + Box(x_box.x(), x_box.y(), x_box.z()), + description.material(x_box.attr<std::string>(_U(material)))); + + // Set envelope volume attributes + envelope_vol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),x_det.visStr()); + + xml_comp_t x_param = x_det.child(_U(param)); + Box box (x_param.x(), x_param.y(), x_param.z()); + Volume box_vol (x_det.nameStr()+"_param", box, description.material(x_param.materialStr())); + xml_dim_t x_trafo = x_param.transformation(); + xml_dim_t x_dim_x = x_trafo.child(_U(dim_x)); + Transform3D start, trafo1; + PlacedVolume pv; + + if ( x_param.hasChild(_U(start)) ) { + start = get_trafo(x_param.child(_U(start))); + } + + trafo1 = get_trafo(x_trafo.child(_U(dim_x))); + + if ( x_trafo.hasChild(_U(dim_y)) && x_trafo.hasChild(_U(dim_z)) ) { + xml_dim_t x_dim_y = x_trafo.child(_U(dim_y)); + xml_dim_t x_dim_z = x_trafo.child(_U(dim_z)); + Transform3D trafo2 = get_trafo(x_dim_y); + Transform3D trafo3 = get_trafo(x_dim_z); + pv = envelope_vol.paramVolume3D(start, box_vol, + x_dim_x.repeat(), trafo1, + x_dim_y.repeat(), trafo2, + x_dim_z.repeat(), trafo3); + } + else if ( x_trafo.hasChild(_U(dim_y)) ) { + xml_dim_t x_dim_y = x_trafo.child(_U(dim_y)); + Transform3D trafo2 = get_trafo(x_dim_y); + pv = envelope_vol.paramVolume2D(start, box_vol, + x_dim_x.repeat(), trafo1, + x_dim_y.repeat(), trafo2); + } + else { + pv = envelope_vol.paramVolume1D(start, box_vol, x_dim_x.repeat(), trafo1); + } + if ( sens.isValid() ) { + sens.setType("calorimeter"); + pv.addPhysVolID("volume", 0); + } + box_vol.setSensitiveDetector(sens); + box_vol.setAttributes(description,x_param.regionStr(),x_param.limitsStr(),x_param.visStr()); + + det.setAttributes(description,envelope_vol,x_det.regionStr(),x_det.limitsStr(),x_det.visStr()); + + // Place the calo inside the world + xml_dim_t x_pos = x_det.position(); + xml_dim_t x_rot = x_det.rotation(); + auto mother = description.pickMotherVolume(det); + Transform3D tr(RotationZYX(x_rot.x(), x_rot.y(), x_rot.z()), + Position(x_pos.x(), x_pos.y(), x_pos.z())); + PlacedVolume envelope_plv = mother.placeVolume(envelope_vol, tr); + envelope_plv.addPhysVolID("system",x_det.id()); + det.setPlacement(envelope_plv); + return det; +} + +DECLARE_DETELEMENT(DD4hep_ParamVolume,create_detector) -- GitLab