From d279562565cb55f97e33f649d91a289f1df51467 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Mon, 20 Nov 2023 17:29:38 +0100 Subject: [PATCH 1/7] add xml for calibration --- ...alBarrel_thetamodulemerged_calibration.xml | 136 ++++++++++++++++++ ..._ECalBarrel_thetamodulemerged_upstream.xml | 136 ++++++++++++++++++ 2 files changed, 272 insertions(+) create mode 100644 FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_calibration.xml create mode 100644 FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_upstream.xml diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_calibration.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_calibration.xml new file mode 100644 index 000000000..be8770709 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_calibration.xml @@ -0,0 +1,136 @@ + + + + + Settings for the inclined EM calorimeter. + The barrel is filled with liquid argon. Passive material includes lead in the middle and steal on the outside, glued together. + Passive plates are inclined by a certain angle from the radial direction. + In between of two passive plates there is a readout. + Space between the plate and readout is of trapezoidal shape and filled with liquid argon. + Definition of sizes, visualization settings, readout and longitudinal segmentation are specified. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 + + + + + + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_upstream.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_upstream.xml new file mode 100644 index 000000000..8c96a5939 --- /dev/null +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_ECalBarrel_thetamodulemerged_upstream.xml @@ -0,0 +1,136 @@ + + + + + Settings for the inclined EM calorimeter. + The barrel is filled with liquid argon. Passive material includes lead in the middle and steal on the outside, glued together. + Passive plates are inclined by a certain angle from the radial direction. + In between of two passive plates there is a readout. + Space between the plate and readout is of trapezoidal shape and filled with liquid argon. + Definition of sizes, visualization settings, readout and longitudinal segmentation are specified. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 + + + + + + system:4,cryo:1,type:3,subtype:3,layer:8,module:11,theta:10 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file From 424db0e484e2319ff68e6e04b1bc38319dd9899b Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Mon, 20 Nov 2023 17:32:44 +0100 Subject: [PATCH 2/7] fix comments inside ALLEGRO xml --- FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml index 49c76c53d..ff4266dd9 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml @@ -3,14 +3,14 @@ xmlns:xs="http://www.w3.org/2001/XMLSchema" xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> - - Master compact file describing the latest developments of the FCCee IDEA detector concept with a LAr calorimeter. + Master compact file describing the latest developments of the FCCee ALLEGRO detector concept. From 849f376169b7a5234be66b146856cd8eba9fde97 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Mon, 20 Nov 2023 17:33:25 +0100 Subject: [PATCH 3/7] fix header ifdef --- .../FCCSWGridModuleThetaMergedHandle_k4geo.h | 4 ++-- .../detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h | 4 ++-- .../detectorSegmentations/FCCSWGridPhiEtaHandle_k4geo.h | 4 ++-- .../include/detectorSegmentations/FCCSWGridPhiEta_k4geo.h | 4 ++-- .../detectorSegmentations/FCCSWGridPhiThetaHandle_k4geo.h | 4 ++-- .../include/detectorSegmentations/FCCSWGridPhiTheta_k4geo.h | 4 ++-- .../include/detectorSegmentations/GridEtaHandle_k4geo.h | 4 ++-- .../include/detectorSegmentations/GridEta_k4geo.h | 4 ++-- .../include/detectorSegmentations/GridThetaHandle_k4geo.h | 4 ++-- .../include/detectorSegmentations/GridTheta_k4geo.h | 4 ++-- 10 files changed, 20 insertions(+), 20 deletions(-) diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h index c0b79ea89..cdbb30223 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMergedHandle_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DD4HEP_DDCORE_GRIDMODULETHETAMERGED_H -#define DD4HEP_DDCORE_GRIDMODULETHETAMERGED_H 1 +#ifndef DETECTORSEGMENTATIONS_FCCSWGRIDMODULETHETAMERGEDHANDLE_K4GEO_H +#define DETECTORSEGMENTATIONS_FCCSWGRIDMODULETHETAMERGEDHANDLE_K4GEO_H // FCCSW #include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h index 796ba400d..98fb122ca 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DETSEGMENTATION_FCCSWGridModuleThetaMerged_k4geo_H -#define DETSEGMENTATION_FCCSWGridModuleThetaMerged_k4geo_H +#ifndef DETECTORSEGMENTATIONS_FCCSWGridModuleThetaMerged_k4geo_H +#define DETECTORSEGMENTATIONS_FCCSWGridModuleThetaMerged_k4geo_H // FCCSW #include "detectorSegmentations/GridTheta_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEtaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEtaHandle_k4geo.h index cd53e01fa..286c97f72 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEtaHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEtaHandle_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DD4HEP_DDCORE_GRIDPHIETA_H -#define DD4HEP_DDCORE_GRIDPHIETA_H 1 +#ifndef DETECTORSEGMENTATIONS_GRIDPHIETAHANDLE_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDPHIETAHANDLE_K4GEO_H // FCCSW #include "detectorSegmentations/FCCSWGridPhiEta_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEta_k4geo.h index 7f80185a3..8c4336262 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiEta_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DETSEGMENTATION_GRIDPHIETA_H -#define DETSEGMENTATION_GRIDPHIETA_H +#ifndef DETECTORSEGMENTATIONS_GRIDPHIETA_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDPHIETA_K4GEO_H // FCCSW #include "detectorSegmentations/GridEta_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiThetaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiThetaHandle_k4geo.h index efbfbf8da..a5546a846 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiThetaHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiThetaHandle_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DD4HEP_DDCORE_GRIDPHITHETA_H -#define DD4HEP_DDCORE_GRIDPHITHETA_H 1 +#ifndef DETECTORSEGMENTATIONS_GRIDPHITHETAHANDLE_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDPHITHETAHANDLE_K4GEO_H // FCCSW #include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiTheta_k4geo.h index cc1d0c31d..aae4bab10 100644 --- a/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/FCCSWGridPhiTheta_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DETSEGMENTATION_GRIDPHITHETA_H -#define DETSEGMENTATION_GRIDPHITHETA_H +#ifndef DETECTORSEGMENTATIONS_GRIDPHITHETA_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDPHITHETA_K4GEO_H // FCCSW #include "detectorSegmentations/GridTheta_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/GridEtaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridEtaHandle_k4geo.h index 961f46025..787c53e23 100644 --- a/detectorSegmentations/include/detectorSegmentations/GridEtaHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/GridEtaHandle_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DD4HEP_DDCORE_GRIDETA_H -#define DD4HEP_DDCORE_GRIDETA_H 1 +#ifndef DETECTORSEGMENTATIONS_GRIDETAHANDLE_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDETAHANDLE_K4GEO_H // FCCSW #include "detectorSegmentations/GridEta_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/GridEta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridEta_k4geo.h index 29cd8dabd..ec4270f79 100644 --- a/detectorSegmentations/include/detectorSegmentations/GridEta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/GridEta_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DETSEGMENTATION_GRIDETA_H -#define DETSEGMENTATION_GRIDETA_H +#ifndef DETECTORSEGMENTATIONS_GRIDETA_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDETA_K4GEO_H #include "DDSegmentation/Segmentation.h" diff --git a/detectorSegmentations/include/detectorSegmentations/GridThetaHandle_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridThetaHandle_k4geo.h index 923406665..9f97d1b63 100644 --- a/detectorSegmentations/include/detectorSegmentations/GridThetaHandle_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/GridThetaHandle_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DD4HEP_DDCORE_GRIDTHETA_H -#define DD4HEP_DDCORE_GRIDTHETA_H 1 +#ifndef DETECTORSEGMENTATIONS_GRIDTHETAHANDLE_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDTHETAHANDLE_K4GEO_H // FCCSW #include "detectorSegmentations/GridTheta_k4geo.h" diff --git a/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h b/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h index 84236c9db..0198e9c7b 100644 --- a/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h +++ b/detectorSegmentations/include/detectorSegmentations/GridTheta_k4geo.h @@ -1,5 +1,5 @@ -#ifndef DETSEGMENTATION_GRIDTHETA_H -#define DETSEGMENTATION_GRIDTHETA_H +#ifndef DETECTORSEGMENTATIONS_GRIDTHETA_K4GEO_H +#define DETECTORSEGMENTATIONS_GRIDTHETA_K4GEO_H #include "DDSegmentation/Segmentation.h" From 52a2b13ed03c0132811d0296ea3132d40625d529 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Mon, 20 Nov 2023 17:34:40 +0100 Subject: [PATCH 4/7] update HCAL readout to theta projectivity for FCCee --- .../FCCee_HCalBarrel_TileCal.xml | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalBarrel_TileCal.xml b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalBarrel_TileCal.xml index 0acecdb1d..714fffbe1 100644 --- a/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalBarrel_TileCal.xml +++ b/FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/FCCee_HCalBarrel_TileCal.xml @@ -45,22 +45,22 @@ - - system:4,layer:5,row:9,eta:9,phi:10 + + system:4,layer:5,row:9,theta:9,phi:10 - - - system:4,layer:5,eta:9,phi:10 + + + system:4,layer:5,theta:9,phi:10 - - - system:4,layer:5,row:9,eta:0,phi:10 + + + system:4,layer:5,row:9,theta:0,phi:10 - + From e4d7edbfd06d8e7d741b5eb7c17b7997f8057dc4 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Mon, 20 Nov 2023 17:37:09 +0100 Subject: [PATCH 5/7] update README for calorimetr --- detector/calorimeter/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/detector/calorimeter/README.md b/detector/calorimeter/README.md index e99fa66f4..f7efc63b6 100644 --- a/detector/calorimeter/README.md +++ b/detector/calorimeter/README.md @@ -12,7 +12,7 @@ The documentation about its usage is [here](../../doc/detector/calorimeter/ECalB Original version taken from [FCCDetectors](https://github.com/HEP-FCC/FCCDetectors/blob/main/Detector/DetFCChhECalInclined/src/ECalBarrelInclined_geo.cpp). ### o1_v02 -New version adapted to the theta segmentation with the possibility to have different number of cell merged per layer. The main difference is that now one has to set `ECalBarrelNumLayers` and `ECalBarrelNumPlanes` in the xml while before it was just dynamically computed based on other xml parameters (also, to avoid silent mistakes, the number from the xml and the one computed dynamically must match). +New version, with module-theta based segmentation. In each layer adjacent cells along theta or module directions can be grouped together, with possibly different merging per layer. This is specified with the `mergedCells_Theta` and `mergedModules` vectors in the `segmentation` tag of the xml file. The baseline grouping in theta is by four in all layers except L1 (the strip layer). The baseline grouping in module direction is by two in all layers. The LAr gap has also been slightly adjusted to bring back the number of modules to 1536 (it was 1545 before). The segmentation class needs to know the number of modules, which is passed via the `nModules` parameter of the `segmentation` tag. To ensure that number of modules and layers (length of the mergedXXX vectors) are consistent with number of modules and layers of the detector, the xml defines `ECalBarrelNumLayers` and `ECalBarrelNumPlanes`, and the c++ file doing the detector construction checks that the number of planes and layers calculated dynamically from other parameters matches that in the xml (if not, the code will crash). ## CaloDisks This sub-detector makes calorimeter endcaps (original and reflected). It is used in ALLEGRO detector concept. From a4be6bcfb3b4aafd4ecfbdb96c4f54960f902e70 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Mon, 20 Nov 2023 17:41:03 +0100 Subject: [PATCH 6/7] add code with helper functions for segmentation --- CMakeLists.txt | 1 + detectorCommon/CMakeLists.txt | 30 + .../include/detectorCommon/DetUtils_k4geo.h | 200 ++++++ detectorCommon/src/DetUtils_k4geo.cpp | 600 ++++++++++++++++++ 4 files changed, 831 insertions(+) create mode 100644 detectorCommon/CMakeLists.txt create mode 100644 detectorCommon/include/detectorCommon/DetUtils_k4geo.h create mode 100644 detectorCommon/src/DetUtils_k4geo.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a13ada8a..384999510 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -54,6 +54,7 @@ find_package( Geant4 REQUIRED ) find_package( LCIO REQUIRED) add_subdirectory(detectorSegmentations) +add_subdirectory(detectorCommon) file(GLOB sources ./detector/tracker/*.cpp diff --git a/detectorCommon/CMakeLists.txt b/detectorCommon/CMakeLists.txt new file mode 100644 index 000000000..4bc3abc9d --- /dev/null +++ b/detectorCommon/CMakeLists.txt @@ -0,0 +1,30 @@ +################################################################################ +# Package: detectorCommon +################################################################################ + +file(GLOB sources src/*.cpp) +add_library(detectorCommon SHARED ${sources}) +target_link_libraries(detectorCommon DD4hep::DDCore DD4hep::DDG4 detectorSegmentations) +target_include_directories(detectorCommon + PUBLIC + $ + $ +) + +file(GLOB headers include/detectorCommon/*.h) +set_target_properties(detectorCommon PROPERTIES PUBLIC_HEADER "${headers}") + +install(TARGETS detectorCommon + EXPORT ${PROJECT_NAME}Targets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_PREFIX}/include/detectorCommon" + COMPONENT dev) + + +dd4hep_generate_rootmap(detectorCommon) + +#include(CTest) +#gaudi_add_test(DumpSimpleBox +# WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} +# FRAMEWORK tests/dumpSimpleBox.py) diff --git a/detectorCommon/include/detectorCommon/DetUtils_k4geo.h b/detectorCommon/include/detectorCommon/DetUtils_k4geo.h new file mode 100644 index 000000000..04c1f93e6 --- /dev/null +++ b/detectorCommon/include/detectorCommon/DetUtils_k4geo.h @@ -0,0 +1,200 @@ +#ifndef DETECTORCOMMON_DETUTILS_H +#define DETECTORCOMMON_DETUTILS_H + +// k4geo +#include "detectorSegmentations/FCCSWGridPhiEta_k4geo.h" +#include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h" +#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h" + +// DD4hep +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Segmentations.h" +#include "DDSegmentation/BitFieldCoder.h" +// Include GridPhiEta from dd4hep once eta calculation is fixed +//#include "DDSegmentation/GridPhiEta.h" +#include "DDSegmentation/CartesianGridXY.h" +#include "DDSegmentation/CartesianGridXYZ.h" +#include "DDSegmentation/PolarGridRPhi.h" + +// Geant +#include "G4Step.hh" + +// CLHEP +#include "CLHEP/Vector/ThreeVector.h" + +#include "TGeoManager.h" + +/** Given a XML element with several daughters with the same name, e.g. + + this method returns the first daughter of type nodeName whose attribute has a given value + e.g. returns when called with (detector, "layer", "name", "1") */ +namespace det { +namespace utils { +dd4hep::xml::Component getNodeByStrAttr(const dd4hep::xml::Handle_t& mother, const std::string& nodeName, + const std::string& attrName, const std::string& attrValue); + +/// try to get attribute with double value, return defaultValue if attribute not found +double getAttrValueWithFallback(const dd4hep::xml::Component& node, const std::string& attrName, + const double& defaultValue); + +/** Retrieves the cellID based on the position of the step and the detector segmentation. + * @param aSeg Handle to the segmentation of the volume. + * @param aStep Step in which particle deposited the energy. + * @param aPreStepPoint Flag indicating if the position of the deposit is the beginning of the step (default) + * or the middle of the step. + */ + +uint64_t cellID(const dd4hep::Segmentation& aSeg, const G4Step& aStep, bool aPreStepPoint = true); + +/** Get number of possible combinations of bit fields for determination of neighbours. + * @param[in] aN number of field names. + * @param[in] aK length of bit fields included for index search. + * return vector of possible combinations of field values. + */ + +std::vector> combinations(int N, int K); + +/** Get number of possible permutations for certain combination of bit field indeces. + * @param[in] aN number of field names. + * return vector of permuations for certain field values. + */ + +std::vector> permutations(int K); + +/** Get true field value of neighbour in cyclic bit field + * @param[in] aCyclicId field value of neighbour + * @param[in] aFieldExtremes Minimal and Maximal values of the fields. + * return field value for neighbour in cyclic bit field + */ + +int cyclicNeighbour(int aCyclicId, std::pair aFieldExtremes); + +/** Get neighbours in many dimensions. + * @param[in] aDecoder Handle to the bitfield decoder. + * @param[in] aFieldNames Names of the fields for which neighbours are found. + * @param[in] aFieldExtremes Minimal and maximal values for the fields. + * @param[in] aCellId ID of cell. + * @param[in] aDiagonal If diagonal neighbours should be included (all combinations of fields). + * return Vector of neighbours. + */ +std::vector neighbours(const dd4hep::DDSegmentation::BitFieldCoder& aDecoder, + const std::vector& aFieldNames, + const std::vector>& aFieldExtremes, + uint64_t aCellId, + const std::vector& aFieldCyclic = {false, false, false, false}, + bool aDiagonal = true); + +/** Special version of the neighbours function for the readout with module and theta merged cells + * Compared to the standard version, it needs a reference to the segmentation class to + * access the number of merged cells per layer. The other parameters and return value are the same + * @param[in] aSeg Reference to the segmentation object. + * @param[in] aDecoder Handle to the bitfield decoder. + * @param[in] aFieldNames Names of the fields for which neighbours are found. + * @param[in] aFieldExtremes Minimal and maximal values for the fields. + * @param[in] aCellId ID of cell. + * @param[in] aDiagonal If diagonal neighbours should be included (all combinations of fields). + * return Vector of neighbours. + */ +std::vector neighbours_ModuleThetaMerged(const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo& aSeg, + const dd4hep::DDSegmentation::BitFieldCoder& aDecoder, + const std::vector& aFieldNames, + const std::vector>& aFieldExtremes, + uint64_t aCellId, + bool aDiagonal = false); + +/** Get minimal and maximal values that can be decoded in the fields of the bitfield. + * @param[in] aDecoder Handle to the bitfield decoder. + * @param[in] aFieldNames Names of the fields for which extremes are found. + * return Vector of pairs (min,max) + */ +std::vector> bitfieldExtremes(const dd4hep::DDSegmentation::BitFieldCoder& aDecoder, + const std::vector& aFieldNames); + +/** Get the half widths of the box envelope (TGeoBBox). + * @param[in] aVolumeId The volume ID. + * return Half-widths of the volume (x,y,z). + */ +CLHEP::Hep3Vector envelopeDimensions(uint64_t aVolumeId); + +/** Get the dimensions of a tube (TGeoConeSeg). + * @param[in] aVolumeId The volume ID. + * return Dimensions of the tube (rmin, rmax, z(half-length)). + */ +CLHEP::Hep3Vector tubeDimensions(uint64_t aVolumeId); + +/** Get the dimensions of a cone (TGeoCone). + * @param[in] aVolumeId The volume ID. + * return Dimensions of the cone (rmin, rmax, z(half-length)). + */ +CLHEP::Hep3Vector coneDimensions(uint64_t aVolumeId); + +/** Get the extrema in pseudorapidity of a tube or cone volume. + * @param[in] aVolumeId The volume ID. + * return Pseudorapidity extrema (eta_min, eta_max). + */ +std::array tubeEtaExtremes(uint64_t aVolumeId); + +/** Get the extrema in pseudorapidity of an envelope. + * @param[in] aVolumeId The volume ID. + * return Pseudorapidity extrema (eta_min, eta_max). + */ +std::array envelopeEtaExtremes(uint64_t aVolumeId); + +/** Get the extrema in pseudorapidity of a volume. First try to match tube or cone, if it fails use an envelope shape. + * @param[in] aVolumeId The volume ID. + * return Pseudorapidity extrema (eta_min, eta_max). + */ +std::array volumeEtaExtremes(uint64_t aVolumeId); + +/** Get the number of cells for the volume and a given Cartesian XY segmentation. + * For an example see: Test/TestReconstruction/tests/options/testcellcountingXYZ.py. + * @warning No offset in segmentation is currently taken into account. + * @param[in] aVolumeId The volume for which the cells are counted. + * @param[in] aSeg Handle to the segmentation of the volume. + * return Array of the number of cells in (X, Y). + */ +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::CartesianGridXY& aSeg); + +/** Get the number of cells for the volume and a given Cartesian XYZ segmentation. + * For an example see: Test/TestReconstruction/tests/options/testcellcountingXYZ.py. + * @warning No offset in segmentation is currently taken into account. + * @param[in] aVolumeId The volume for which the cells are counted. + * @param[in] aSeg Handle to the segmentation of the volume. + * return Array of the number of cells in (X, Y, Z). + */ +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::CartesianGridXYZ& aSeg); + +/** Get the number of cells for the volume and a given Phi-Eta / Phi-Theta / Module-Theta segmentation. + * It is assumed that the volume has a cylindrical shape (and full azimuthal coverage) + * and that it is centred at (0,0,0). + * For an example see: Test/TestReconstruction/tests/options/testcellcountingPhiEta.py. + * @warning No offset in segmentation is currently taken into account. + * @param[in] aVolumeId The volume for which the cells are counted. + * @param[in] aSeg Handle to the segmentation of the volume. + * return Array of the number of cells in (phi, eta) / (phi, theta) / (module, theta) and the minimum eta / theta ID. + */ +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridPhiEta_k4geo& aSeg); +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo& aSeg); +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo& aSeg); + +/** Get the number of cells for the volume and a given R-phi segmentation. + * It is assumed that the volume has a cylindrical shape - TGeoTube (and full azimuthal coverage) + * and that it is centred at (0,0,0). + * For an example see: Test/TestReconstruction/tests/options/testcellcountingRPhi.py. + * @warning No offset in segmentation is currently taken into account. + * @param[in] aVolumeId The volume for which the cells are counted. + * @param[in] aSeg Handle to the segmentation of the volume. + * return Array of the number of cells in (r, phi). + */ +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::PolarGridRPhi& aSeg); + +/** Get the number of the volumes containing a given name. + * For an example see: Test/TestReconstruction/tests/options/testcellcountingXYZ.py. + * @param[in] aHighestVolume The top volume in the geometry. + * @param[in] aMatchName Name (or its part) of the volume. + * return Number of the volumes. + */ +unsigned int countPlacedVolumes(TGeoVolume* aHighestVolume, const std::string& aMatchName); +} +} +#endif /* DETCOMMON_DETUTILS_H */ diff --git a/detectorCommon/src/DetUtils_k4geo.cpp b/detectorCommon/src/DetUtils_k4geo.cpp new file mode 100644 index 000000000..711118248 --- /dev/null +++ b/detectorCommon/src/DetUtils_k4geo.cpp @@ -0,0 +1,600 @@ +#include "detectorCommon/DetUtils_k4geo.h" + +// DD4hep +#include "DDG4/Geant4Mapping.h" +#include "DDG4/Geant4VolumeManager.h" + +// Geant +#include "G4NavigationHistory.hh" + +// ROOT +#include "TGeoBBox.h" + +#ifdef HAVE_GEANT4_UNITS +#define MM_2_CM 1.0 +#else +#define MM_2_CM 0.1 +#endif + +#include + +#include + +namespace det { +namespace utils { +dd4hep::xml::Component getNodeByStrAttr(const dd4hep::xml::Handle_t& mother, const std::string& nodeName, + const std::string& attrName, const std::string& attrValue) { + for (dd4hep::xml::Collection_t xCompColl(mother, nodeName.c_str()); nullptr != xCompColl; ++xCompColl) { + if (xCompColl.attr(attrName.c_str()) == attrValue) { + return static_cast(xCompColl); + } + } + // in case there was no xml daughter with matching name + return dd4hep::xml::Component(nullptr); +} + +double getAttrValueWithFallback(const dd4hep::xml::Component& node, const std::string& attrName, + const double& defaultValue) { + if (node.hasAttr(_Unicode(attrName.c_str()))) { + return node.attr(attrName.c_str()); + } else { + return defaultValue; + } +} + +uint64_t cellID(const dd4hep::Segmentation& aSeg, const G4Step& aStep, bool aPreStepPoint) { + dd4hep::sim::Geant4VolumeManager volMgr = dd4hep::sim::Geant4Mapping::instance().volumeManager(); + dd4hep::VolumeID volID = volMgr.volumeID(aStep.GetPreStepPoint()->GetTouchable()); + if (aSeg.isValid()) { + G4ThreeVector global; + if (aPreStepPoint) { + global = aStep.GetPreStepPoint()->GetPosition(); + } else { + global = 0.5 * (aStep.GetPreStepPoint()->GetPosition() + aStep.GetPostStepPoint()->GetPosition()); + } + G4ThreeVector local = + aStep.GetPreStepPoint()->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(global); + dd4hep::Position loc(local.x() * MM_2_CM, local.y() * MM_2_CM, local.z() * MM_2_CM); + dd4hep::Position glob(global.x() * MM_2_CM, global.y() * MM_2_CM, global.z() * MM_2_CM); + dd4hep::VolumeID cID = aSeg.cellID(loc, glob, volID); + return cID; + } + return volID; +} + +std::vector> combinations(int N, int K) { + std::vector> indexes; + std::string bitmask(K, 1); // K leading 1's + bitmask.resize(N, 0); // N-K trailing 0's + // permute bitmask + do { + std::vector tmp; + for (int i = 0; i < N; ++i) { // [0..N-1] integers + if (bitmask[i]) { + tmp.push_back(i); + } + } + indexes.push_back(tmp); + } while (std::prev_permutation(bitmask.begin(), bitmask.end())); + return indexes; +} + +std::vector> permutations(int K) { + std::vector> indexes; + int N = pow(2, K); // number of permutations with repetition of 2 numbers (-1,1) + for (int i = 0; i < N; i++) { + // permutation = binary representation of i + std::vector tmp; + tmp.assign(K, 0); + uint res = i; + // dec -> bin + for (int j = 0; j < K; j++) { + tmp[K - 1 - j] = -1 + 2 * (res % 2); + res = floor(res / 2); + } + indexes.push_back(tmp); + } + return indexes; +} + +// use it for field module/phi +int cyclicNeighbour(int aCyclicId, std::pair aFieldExtremes) { + int nBins = aFieldExtremes.second - aFieldExtremes.first + 1; + if (aCyclicId < aFieldExtremes.first) { + return (aFieldExtremes.second + 1 - ((aFieldExtremes.first - aCyclicId) % nBins) ); + } else if (aCyclicId > aFieldExtremes.second) { + return ( ((aCyclicId - aFieldExtremes.first) % nBins) + aFieldExtremes.first) ; + } + return aCyclicId; +} + +std::vector neighbours(const dd4hep::DDSegmentation::BitFieldCoder& aDecoder, + const std::vector& aFieldNames, + const std::vector>& aFieldExtremes, uint64_t aCellId, + const std::vector& aFieldCyclic, bool aDiagonal) { + std::vector neighbours; + dd4hep::DDSegmentation::CellID cID = aCellId; + for (uint itField = 0; itField < aFieldNames.size(); itField++) { + const auto& field = aFieldNames[itField]; + // note: get(..) returns a FieldID i.e. a int64_t + // while CellID (used previously) is a uint64_t + // similarly, the second argument of set(..) + // is signed (FieldID) rather than unsigned (CellID) + int id = aDecoder.get(cID,field); + if (aFieldCyclic[itField]) { + aDecoder[field].set(cID, cyclicNeighbour(id - 1, aFieldExtremes[itField])); + neighbours.emplace_back(cID); + aDecoder[field].set(cID, cyclicNeighbour(id + 1, aFieldExtremes[itField])); + neighbours.emplace_back(cID); + } else { + if (id > aFieldExtremes[itField].first) { + aDecoder[field].set(cID, id - 1); + neighbours.emplace_back(cID); + } + if (id < aFieldExtremes[itField].second) { + aDecoder[field].set(cID, id + 1); + neighbours.emplace_back(cID); + } + } + aDecoder[field].set(cID, id); + } + if (aDiagonal) { + std::vector fieldIds; // initial IDs + fieldIds.assign(aFieldNames.size(), 0); + // for each field get current Id + for (uint iField = 0; iField < aFieldNames.size(); iField++) { + const auto& field = aFieldNames[iField]; + fieldIds[iField] = aDecoder.get(cID, field); + } + for (uint iLength = aFieldNames.size(); iLength > 1; iLength--) { + // get all combinations for a given length + const auto& indexes = combinations(aFieldNames.size(), iLength); + for (uint iComb = 0; iComb < indexes.size(); iComb++) { + // for current combination get all permutations of +- 1 operation on IDs + const auto& calculation = permutations(iLength); + // do the increase/decrease of bitfield + for (uint iCalc = 0; iCalc < calculation.size(); iCalc++) { + // set new Ids for each field combination + bool add = true; + for (uint iField = 0; iField < indexes[iComb].size(); iField++) { + if (aFieldCyclic[indexes[iComb][iField]]) { + aDecoder[aFieldNames[indexes[iComb][iField]]].set(cID, cyclicNeighbour(fieldIds[indexes[iComb][iField]] + calculation[iCalc][iField], + aFieldExtremes[indexes[iComb][iField]]) ); + } else if ((calculation[iCalc][iField] > 0 && + fieldIds[indexes[iComb][iField]] < aFieldExtremes[indexes[iComb][iField]].second) || + (calculation[iCalc][iField] < 0 && + fieldIds[indexes[iComb][iField]] > aFieldExtremes[indexes[iComb][iField]].first)) { + aDecoder[aFieldNames[indexes[iComb][iField]]].set(cID, fieldIds[indexes[iComb][iField]] + calculation[iCalc][iField]); + } else { + add = false; + } + } + // add new cellId to neighbours (unless it's beyond extrema) + if (add) { + neighbours.emplace_back(cID); + } + // reset ids + for (uint iField = 0; iField < indexes[iComb].size(); iField++) { + aDecoder[aFieldNames[indexes[iComb][iField]]].set(cID, fieldIds[indexes[iComb][iField]]); + } + } + } + } + } + return neighbours; +} + +// use it for module-theta merged readout (FCCSWGridModuleThetaMerged_k4geo) +std::vector neighbours_ModuleThetaMerged(const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo& aSeg, + const dd4hep::DDSegmentation::BitFieldCoder& aDecoder, + const std::vector& aFieldNames, + const std::vector>& aFieldExtremes, + uint64_t aCellId, + bool aDiagonal) { + + std::vector neighbours; + + // check that field names and extremes have the proper length + if (aFieldNames.size() != 3 || aFieldExtremes.size()!=3) { + std::cout << "ERROR: the vectors aFieldNames and aFieldSizes should be of length = 3, corresponding to the theta/module/layer fields" << std::endl; + std::cout << "ERROR: will return empty neighbour map" << std::endl; + return neighbours; + } + + // find index of layer, module and theta in field extremes vector + int idModuleField(-1); + int idThetaField(-1); + int idLayerField(-1); + for (uint itField = 0; itField < aFieldNames.size(); itField++) { + if (aFieldNames[itField] == aSeg.fieldNameModule()) + idModuleField = (int) itField; + else if (aFieldNames[itField] == aSeg.fieldNameTheta()) + idThetaField = (int) itField; + else if (aFieldNames[itField] == aSeg.fieldNameLayer()) + idLayerField = (int) itField; + + } + if (idModuleField < 0) { + std::cout << "WARNING: module field " << aSeg.fieldNameModule() << " not found in aFieldNames vector" << std::endl; + std::cout << "WARNING: will return empty neighbour map" << std::endl; + return neighbours; + } + if (idThetaField < 0) { + std::cout << "WARNING: theta field " << aSeg.fieldNameTheta() << " not found in aFieldNames vector" << std::endl; + std::cout << "WARNING: will return empty neighbour map" << std::endl; + return neighbours; + } + if (idLayerField < 0) { + std::cout << "WARNING: layer field " << aSeg.fieldNameLayer() << " not found in aFieldNames vector" << std::endl; + std::cout << "WARNING: will return empty neighbour map" << std::endl; + return neighbours; + } + + // retrieve layer/module/theta of cell under study + int layer_id = aDecoder.get(aCellId, aFieldNames[idLayerField]); + int module_id = aDecoder.get(aCellId, aFieldNames[idModuleField]); + int theta_id = aDecoder.get(aCellId, aFieldNames[idThetaField]); + + // now find the neighbours + dd4hep::DDSegmentation::CellID cID = aCellId; + + // for neighbours across different layers, we have to take into + // account that merging along module and/or theta could be different + // so one cell in layer N could be neighbour to several in layer N+-1 + // The cells are classified in a different way whether they are + // direct neighbours (common surface), diagonal neighbours (common edge or vertex) + // or neither. + // To decide this, we need to check how the cells are related in both directions: + // neighbours (edge at least partially in common), diagonal neigbours (common vertex), + // none + int neighbourTypeModuleDir; // 0: not neighbour; 1: diagonal neighbour; 2: neighbour in module direction + int neighbourTypeThetaDir; // 0: not neighbour; 1: diagonal neighbour; 2: neighbour in module direction + + for (int deltaLayer = -1; deltaLayer<2; deltaLayer+=2) { + + // no neighbours in layer N-1 for innermost layer + if (layer_id == aFieldExtremes[idLayerField].first && deltaLayer<0) continue; + // and in layer N+1 for outermost layer + if (layer_id == aFieldExtremes[idLayerField].second && deltaLayer>0) continue; + + // set layer field of neighbour cell + aDecoder.set(cID, aSeg.fieldNameLayer(), layer_id + deltaLayer); + + // find the neighbour(s) in module and theta + // if the numbers of module (theta) merged cells across 2 layers are the + // same then we just take the same module (theta) ID + // otherwise, we need to do some math to account for the different mergings + // note: even if number of merged cells in layer-1 is larger, a cell + // in layer could neighbour more than one cell in layer-1 if the merged + // cells are not aligned, for example if cells are grouped by 3 in a layer + // and by 4 in the next one, cell 435 in the former (which groups together + // 435-436-437) will be neighbour to cells 432 and 436 of the latter + // this might introduce duplicates, we will remove them later + // another issue is that it could add spurious cells beyond the maximum module number + // to prevent this we would need to know the max module number in layer -1 + // which would require modifying this function passing the extrema for all layers + // instead of the extrema only for a certain layer + // this border effect is also present in the original method.. + for (int i=-1; i <= aSeg.mergedThetaCells(layer_id); i++) { + int theta_id_neighbour = (theta_id + i) - ((theta_id + i) % aSeg.mergedThetaCells(layer_id+deltaLayer)); + if (theta_id_neighbour >= theta_id && theta_id_neighbour < (theta_id + aSeg.mergedThetaCells(layer_id))) neighbourTypeThetaDir = 2; + else if (theta_id_neighbour < theta_id && theta_id_neighbour > (theta_id - aSeg.mergedThetaCells(layer_id+deltaLayer))) neighbourTypeThetaDir = 2; + else if (theta_id_neighbour == (theta_id + aSeg.mergedThetaCells(layer_id))) neighbourTypeThetaDir = 1; + else if (theta_id_neighbour == (theta_id - aSeg.mergedThetaCells(layer_id+deltaLayer))) neighbourTypeThetaDir = 1; + else neighbourTypeThetaDir = 0; + + // if there is no point of contact along theta, no need to check also for module direction + if (neighbourTypeThetaDir == 0) continue; + // if we are not considering diagonal neighbours, and cells in theta have only an edge in common, then skip + if (!aDiagonal && neighbourTypeThetaDir == 1) continue; + // otherwise, check status along module direction + for (int j=-1; j <= aSeg.mergedModules(layer_id); j++) { + int module_id_neighbour = (module_id + j) - ((module_id + j) % aSeg.mergedModules(layer_id+deltaLayer)); + int module_id_neighbour_cyclic = cyclicNeighbour(module_id_neighbour, + aFieldExtremes[idModuleField] + ); + + if (module_id_neighbour >= module_id && module_id_neighbour < (module_id + aSeg.mergedModules(layer_id))) neighbourTypeModuleDir = 2; + else if (module_id_neighbour < module_id && module_id_neighbour > (module_id - aSeg.mergedModules(layer_id+deltaLayer))) neighbourTypeModuleDir = 2; + else if (module_id_neighbour == (module_id + aSeg.mergedModules(layer_id))) neighbourTypeModuleDir = 1; + else if (module_id_neighbour == (module_id - aSeg.mergedModules(layer_id+deltaLayer))) neighbourTypeModuleDir = 1; + else neighbourTypeModuleDir = 0; + + // if there is no point of contact along module, then skip + if (neighbourTypeModuleDir == 0) continue; + // otherwise: if neighbours along both module and theta, or along one of the two + // and we also consider diagonal neighbours, then add cells to list of neighbours + if ( (neighbourTypeModuleDir == 2 && neighbourTypeThetaDir==2) || + (aDiagonal && neighbourTypeThetaDir > 0 && neighbourTypeModuleDir>0) ) { + aDecoder.set(cID, aSeg.fieldNameModule(), module_id_neighbour_cyclic); + aDecoder.set(cID, aSeg.fieldNameTheta(), theta_id_neighbour); + neighbours.emplace_back(cID); + } + } + } + } + // reset cellID + // aDecoder.set(cID, aSeg.fieldNameModule(), module_id); + // aDecoder.set(cID, aSeg.fieldNameTheta(), theta_id); + + // for neighbours in module/theta direction at same layer_id, do +-nMergedCells instead of +-1 + aDecoder.set(cID, aSeg.fieldNameLayer(), layer_id); + // loop over theta cells + for (int i=-1; i<=1; i++) { + // calculate theta_id of neighbour + int theta_id_neighbour = theta_id + i*aSeg.mergedThetaCells(layer_id); + // check that it is within the ranges + if ( + (theta_id_neighbour < aFieldExtremes[idThetaField].first) || + (theta_id_neighbour > aFieldExtremes[idThetaField].second) + ) continue; + // set theta_id of cell ID + aDecoder[aSeg.fieldNameTheta()].set(cID, theta_id + i*aSeg.mergedThetaCells(layer_id)); + // loop over modules + for (int j=-1; j<=1; j++) { + // skip the cell under study (i==j==0) + if (i==0 && j==0) continue; + // calculate module_id of neighbour + int newid = cyclicNeighbour(module_id + j*aSeg.mergedModules(layer_id), aFieldExtremes[idModuleField]); + newid -= (newid % aSeg.mergedModules(layer_id)); + // set module_if of cell ID + aDecoder[aSeg.fieldNameModule()].set(cID, newid); + // add cell to neighbour list + if ( i==0 || j==0 || aDiagonal ) { // first two conditions correspond to non diagonal neighbours + neighbours.emplace_back(cID); + } + } + } + + // remove duplicates + std::unordered_set s; + for (uint64_t i : neighbours) + s.insert(i); + neighbours.assign( s.begin(), s.end() ); + return neighbours; +} + +std::vector> bitfieldExtremes(const dd4hep::DDSegmentation::BitFieldCoder& aDecoder, + const std::vector& aFieldNames) { + std::vector> extremes; + int width = 0; + for (const auto& field : aFieldNames) { + width = aDecoder[field].width(); + if (aDecoder[field].isSigned()) { + extremes.emplace_back(std::make_pair(-(1 << (width - 1)), (1 << (width - 1)) - 1)); + } else { + extremes.emplace_back(std::make_pair(0, (1 << width) - 1)); + } + } + return extremes; +} + +CLHEP::Hep3Vector envelopeDimensions(uint64_t aVolumeId) { + dd4hep::VolumeManager volMgr = dd4hep::Detector::getInstance().volumeManager(); + auto pvol = volMgr.lookupVolumePlacement(aVolumeId); + auto solid = pvol.volume().solid(); + // get the envelope of the shape + TGeoBBox* box = (dynamic_cast(solid.ptr())); + // get half-widths + return CLHEP::Hep3Vector(box->GetDX(), box->GetDY(), box->GetDZ()); +} + +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::CartesianGridXY& aSeg) { + // // get half-widths + auto halfSizes = envelopeDimensions(aVolumeId); + // get segmentation cell widths + double xCellSize = aSeg.gridSizeX(); + double yCellSize = aSeg.gridSizeY(); + // calculate number of cells, the middle cell is centred at 0 (no offset) + uint cellsX = ceil((halfSizes.x() - xCellSize / 2.) / xCellSize) * 2 + 1; + uint cellsY = ceil((halfSizes.y() - yCellSize / 2.) / yCellSize) * 2 + 1; + return {cellsX, cellsY}; +} + +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::CartesianGridXYZ& aSeg) { + // // get half-widths + auto halfSizes = envelopeDimensions(aVolumeId); + // get segmentation cell widths + double xCellSize = aSeg.gridSizeX(); + double yCellSize = aSeg.gridSizeY(); + double zCellSize = aSeg.gridSizeZ(); + // calculate number of cells, the middle cell is centred at 0 (no offset) + uint cellsX = ceil((halfSizes.x() - xCellSize / 2.) / xCellSize) * 2 + 1; + uint cellsY = ceil((halfSizes.y() - yCellSize / 2.) / yCellSize) * 2 + 1; + uint cellsZ = ceil((halfSizes.z() - zCellSize / 2.) / zCellSize) * 2 + 1; + return {cellsX, cellsY, cellsZ}; +} + +CLHEP::Hep3Vector tubeDimensions(uint64_t aVolumeId) { + dd4hep::VolumeManager volMgr = dd4hep::Detector::getInstance().volumeManager(); + auto pvol = volMgr.lookupVolumePlacement(aVolumeId); + auto solid = pvol.volume().solid(); + + // get the envelope of the shape + TGeoTubeSeg* tube = (dynamic_cast(solid.ptr())); + if (tube == nullptr) { + return CLHEP::Hep3Vector(0, 0, 0); + } + // get half-widths + return CLHEP::Hep3Vector(tube->GetRmin(), tube->GetRmax(), tube->GetDZ()); +} + +CLHEP::Hep3Vector coneDimensions(uint64_t aVolumeId) { + dd4hep::VolumeManager volMgr = dd4hep::Detector::getInstance().volumeManager(); + auto pvol = volMgr.lookupVolumePlacement(aVolumeId); + auto solid = pvol.volume().solid(); + // get the envelope of the shape + TGeoCone* cone = (dynamic_cast(solid.ptr())); + if (cone == nullptr) { + return CLHEP::Hep3Vector(0, 0, 0); + } + // get half-widths + return CLHEP::Hep3Vector(cone->GetRmin1(), cone->GetRmax1(), cone->GetDZ()); +} + +std::array tubeEtaExtremes(uint64_t aVolumeId) { + auto sizes = tubeDimensions(aVolumeId); + if (sizes.mag() == 0) { + // if it is not a cylinder maybe it is a cone (same calculation for extremes) + sizes = coneDimensions(aVolumeId); + if (sizes.mag() == 0) { + return {0, 0}; + } + } + // eta segmentation calculate maximum eta from the inner radius (no offset is taken into account) + double maxEta = 0; + double minEta = 0; + + double rIn = sizes.x(); + double rOut = sizes.y(); + double dZ = sizes.z(); + + // check if it is a cylinder centred at z=0 + dd4hep::VolumeManager volMgr = dd4hep::Detector::getInstance().volumeManager(); + auto detelement = volMgr.lookupDetElement(aVolumeId); + const auto& transformMatrix = detelement.nominal().worldTransformation(); + double outGlobal[3]; + double inLocal[] = {0, 0, 0}; // to get middle of the volume + transformMatrix.LocalToMaster(inLocal, outGlobal); + double zCenter = outGlobal[2]; + if (fabs(zCenter) < 1e-10) { + // this assumes cylinder centred at z=0 + maxEta = CLHEP::Hep3Vector(rIn, 0, dZ).eta(); + minEta = -maxEta; + } else { + maxEta = std::max( + CLHEP::Hep3Vector(rIn, 0, zCenter+dZ).eta(), + CLHEP::Hep3Vector(rOut, 0, zCenter+dZ).eta() + ); + minEta = std::min( + CLHEP::Hep3Vector(rIn, 0, zCenter-dZ).eta(), + CLHEP::Hep3Vector(rOut, 0, zCenter-dZ).eta() + ); + } + return {minEta, maxEta}; +} + +std::array envelopeEtaExtremes (uint64_t aVolumeId) { + dd4hep::VolumeManager volMgr = dd4hep::Detector::getInstance().volumeManager(); + auto detelement = volMgr.lookupDetElement(aVolumeId); + const auto& transformMatrix = detelement.nominal().worldTransformation(); + // calculate values of eta in all possible corners of the envelope + auto dim = envelopeDimensions(aVolumeId); + double minEta = 0; + double maxEta = 0; + for (uint i = 0; i < 8; i++) { + // coefficients to get all combinations of corners + int iX = -1 + 2 * ((i / 2) % 2); + int iY = -1 + 2 * (i % 2); + int iZ = -1 + 2 * (i / 4); + double outDimGlobal[3]; + double inDimLocal[] = {iX * dim.x(), iY * dim.y(), iZ * dim.z()}; + transformMatrix.LocalToMaster(inDimLocal, outDimGlobal); + double eta = CLHEP::Hep3Vector(outDimGlobal[0], outDimGlobal[1], outDimGlobal[2]).eta(); + if (i == 0) { + minEta = eta; + maxEta = eta; + } + if (eta < minEta) { + minEta = eta; + } + if (eta > maxEta) { + maxEta = eta; + } + } + return {minEta, maxEta}; +} + +std::array volumeEtaExtremes(uint64_t aVolumeId) { + // try if volume is a cylinder/disc + auto etaExtremes = tubeEtaExtremes(aVolumeId); + if (etaExtremes[0] != 0 or etaExtremes[1] != 0) { + return etaExtremes; + } else { + return envelopeEtaExtremes(aVolumeId); + } +} + +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridPhiEta_k4geo& aSeg) { + // get segmentation number of bins in phi + uint phiCellNumber = aSeg.phiBins(); + // get segmentation cell width in eta + double etaCellSize = aSeg.gridSizeEta(); + // get min and max eta of the volume + auto etaExtremes = volumeEtaExtremes(aVolumeId); + // calculate the number of eta volumes + // max - min = full eta range, - size = not counting the middle cell centred at 0, + 1 to account for that cell + uint cellsEta = ceil(( etaExtremes[1] - etaExtremes[0] - etaCellSize ) / 2 / etaCellSize) * 2 + 1; + uint minEtaID = int(floor((etaExtremes[0] + 0.5 * etaCellSize - aSeg.offsetEta()) / etaCellSize)); + return {phiCellNumber, cellsEta, minEtaID}; +} + +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridPhiTheta_k4geo& aSeg) { + uint phiCellNumber = aSeg.phiBins(); + double thetaCellSize = aSeg.gridSizeTheta(); + double thetaOffset = aSeg.offsetTheta(); + auto etaExtremes = volumeEtaExtremes(aVolumeId); + double thetaMin = 2.*atan(exp(-etaExtremes[1])); + double thetaMax = 2.*atan(exp(-etaExtremes[0])); + // debug + std::cout << "volumeID = " << aVolumeId << std::endl; + std::cout << "thetaMin, thetaMax = " << thetaMin << " " << thetaMax << std::endl; + // + uint minThetaID = int(floor((thetaMin + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); + uint maxThetaID = int(floor((thetaMax + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); + uint nThetaCells = 1 + maxThetaID - minThetaID; + return {phiCellNumber, nThetaCells, minThetaID}; +} + +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo& aSeg) { + + const dd4hep::DDSegmentation::BitFieldCoder* aDecoder = aSeg.decoder(); + int nLayer = aDecoder->get(aVolumeId, aSeg.fieldNameLayer()); + // + 0.5 to avoid integer division + uint nModules = aSeg.nModules() / aSeg.mergedModules(nLayer); + if (aSeg.nModules() % aSeg.mergedModules(nLayer) != 0) nModules++; + // get minimum and maximum theta of volume + auto etaExtremes = volumeEtaExtremes(aVolumeId); + double thetaMin = 2.*atan(exp(-etaExtremes[1])); + double thetaMax = 2.*atan(exp(-etaExtremes[0])); + + // convert to minimum and maximum theta bins + double thetaCellSize = aSeg.gridSizeTheta(); + double thetaOffset = aSeg.offsetTheta(); + uint minThetaID = int(floor((thetaMin + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); + uint maxThetaID = int(floor((thetaMax + 0.5 * thetaCellSize - thetaOffset) / thetaCellSize)); + // correct minThetaID and maxThetaID for merging + uint mergedThetaCells = aSeg.mergedThetaCells(nLayer); + minThetaID -= (minThetaID % mergedThetaCells); + maxThetaID -= (maxThetaID % mergedThetaCells); + uint nThetaCells = 1 + (maxThetaID - minThetaID)/ mergedThetaCells; + return {nModules, nThetaCells, minThetaID}; +} + +std::array numberOfCells(uint64_t aVolumeId, const dd4hep::DDSegmentation::PolarGridRPhi& aSeg) { + // get half-widths, + auto tubeSizes = tubeDimensions(aVolumeId); + // get segmentation cell width + double rCellSize = aSeg.gridSizeR(); + double phiCellSize = aSeg.gridSizePhi(); + uint cellsRout = ceil(tubeSizes.y() / rCellSize); + uint cellsRin = floor(tubeSizes.x() / rCellSize); + uint cellsR = cellsRout - cellsRin; + uint cellsPhi = ceil(2 * M_PI / phiCellSize); + return {cellsR, cellsPhi}; +} + +unsigned int countPlacedVolumes(TGeoVolume* aHighestVolume, const std::string& aMatchName) { + int numberOfPlacedVolumes = 0; + TGeoNode* node; + TGeoIterator next(aHighestVolume); + while ((node = next())) { + std::string currentNodeName = node->GetName(); + if (currentNodeName.find(aMatchName) != std::string::npos) { + ++numberOfPlacedVolumes; + } + } + return numberOfPlacedVolumes; +} +} +} From 0fd34afc1d26a92dbfe382cf113db21a7eae8805 Mon Sep 17 00:00:00 2001 From: Giovanni Marchiori Date: Mon, 11 Dec 2023 12:34:14 +0100 Subject: [PATCH 7/7] implement suggestions by Andre Pleier --- detectorCommon/CMakeLists.txt | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/detectorCommon/CMakeLists.txt b/detectorCommon/CMakeLists.txt index 4bc3abc9d..0982e0a4c 100644 --- a/detectorCommon/CMakeLists.txt +++ b/detectorCommon/CMakeLists.txt @@ -3,11 +3,11 @@ ################################################################################ file(GLOB sources src/*.cpp) -add_library(detectorCommon SHARED ${sources}) +add_dd4hep_plugin(detectorCommon SHARED ${sources}) target_link_libraries(detectorCommon DD4hep::DDCore DD4hep::DDG4 detectorSegmentations) target_include_directories(detectorCommon PUBLIC - $ + $ $ ) @@ -18,13 +18,4 @@ install(TARGETS detectorCommon EXPORT ${PROJECT_NAME}Targets RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib - PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_PREFIX}/include/detectorCommon" - COMPONENT dev) - - -dd4hep_generate_rootmap(detectorCommon) - -#include(CTest) -#gaudi_add_test(DumpSimpleBox -# WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} -# FRAMEWORK tests/dumpSimpleBox.py) + PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/detectorCommon" COMPONENT dev)