Merge branch 'key4hep:master' into DCH_digi

CMakeLists.txt

@@ -5,7 +5,7 @@ project(${PackageName})
 
 # please keep this layout for version setting, used by the automatic tagging script
 set(${PackageName}_VERSION_MAJOR 0)
-set(${PackageName}_VERSION_MINOR 2)
+set(${PackageName}_VERSION_MINOR 3)
 set(${PackageName}_VERSION_PATCH 0)
 
 set(${PackageName}_VERSION "${PACKAGE_VERSION_MAJOR}.${PACKAGE_VERSION_MINOR}.${PACKAGE_VERSION_PATCH}")
@@ -24,35 +24,9 @@ find_package(Gaudi)
 find_package(GSL)
 #---------------------------------------------------------------
 
+include(cmake/Key4hepConfig.cmake)
 include(GNUInstallDirs)
 
-# Set up C++ Standard
-# ``-DCMAKE_CXX_STANDARD=<standard>`` when invoking CMake
-set(CMAKE_CXX_STANDARD 17 CACHE STRING "")
-
-if(NOT CMAKE_CXX_STANDARD MATCHES "14|17|20")
-  message(FATAL_ERROR "Unsupported C++ standard: ${CMAKE_CXX_STANDARD}")
-endif()
-
-# use, i.e. don't skip the full RPATH for the build tree
-set(CMAKE_SKIP_BUILD_RPATH FALSE)
-
-# when building, don't use the install RPATH already
-# (but later on when installing)
-set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
-
-set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
-
-# add the automatically determined parts of the RPATH
-# which point to directories outside the build tree to the install RPATH
-set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
-
-# the RPATH to be used when installing, but only if it's not a system directory
-list(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir)
-if("${isSystemDir}" STREQUAL "-1")
-    set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
-endif("${isSystemDir}" STREQUAL "-1")
-
 include(CTest)
 function(set_test_env _testname)
 	set_property(TEST ${test_name} APPEND PROPERTY ENVIRONMENT "ROOT_INCLUDE_PATH=$<$<TARGET_EXISTS:podio::podio>:$<TARGET_FILE_DIR:podio::podio>/../include>:$<$<TARGET_EXISTS:EDM4HEP::edm4hep>:$<TARGET_FILE_DIR:EDM4HEP::edm4hep>/../include>:$ENV{ROOT_INCLUDE_PATH}")

Tracking/components/PlotTrackHitResiduals.cpp

@@ -6,7 +6,7 @@
 // edm4hep
 #include "edm4hep/MCParticleCollection.h"
 #include "edm4hep/TrackCollection.h"
-#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/TrackMCParticleLinkCollection.h"
 #include "edm4hep/SimTrackerHitCollection.h"
 
 // marlin
@@ -18,6 +18,17 @@
 // k4FWCore
 #include "k4FWCore/Consumer.h"
 
+#include "GAUDI_VERSION.h"
+
+#if GAUDI_MAJOR_VERSION < 39
+namespace Gaudi::Accumulators {
+  template <unsigned int ND, atomicity Atomicity = atomicity::full, typename Arithmetic = double>
+  using StaticHistogram =
+      Gaudi::Accumulators::HistogramingCounterBase<ND, Atomicity, Arithmetic, naming::histogramString,
+                                                   HistogramingAccumulator>;
+}
+#endif
+
 #include <string>
 
 /** @class PlotTrackHitDistances
@@ -29,7 +40,7 @@
  */
 
 struct PlotTrackHitDistances final
-  : k4FWCore::Consumer<void(const edm4hep::SimTrackerHitCollection&, const edm4hep::MCRecoTrackParticleAssociationCollection&)> {
+  : k4FWCore::Consumer<void(const edm4hep::SimTrackerHitCollection&, const edm4hep::TrackMCParticleLinkCollection&)> {
   PlotTrackHitDistances(const std::string& name, ISvcLocator* svcLoc)
       : Consumer(
             name, svcLoc,
@@ -38,11 +49,11 @@ struct PlotTrackHitDistances final
             KeyValues("InputTracksFromGenParticlesAssociation", {"TracksFromGenParticlesAssociation"}),
             }) {}
 
-  void operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits, const edm4hep::MCRecoTrackParticleAssociationCollection& trackParticleAssociations) const override {
+  void operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits, const edm4hep::TrackMCParticleLinkCollection& trackParticleAssociations) const override {
 
     for (const auto& trackParticleAssociation : trackParticleAssociations) {
-      auto genParticle = trackParticleAssociation.getSim();
-      auto track = trackParticleAssociation.getRec();
+      auto genParticle = trackParticleAssociation.getTo();
+      auto track = trackParticleAssociation.getFrom();
       edm4hep::TrackState trackStateAtIP;
       bool found_trackStateAtIP = false;
       for (const auto& trackState : track.getTrackStates()) {
@@ -74,7 +85,7 @@ struct PlotTrackHitDistances final
     return;
   }
   Gaudi::Property<float> m_Bz{this, "Bz", 2., "Z component of the (assumed constant) magnetic field in Tesla."};
-  mutable Gaudi::Accumulators::Histogram<1> m_residualHist{this, "track_hits_distance_closest_approach", "Track-hit Distances", {100, 0, 1, "Distance [mm];Entries"}};
+  mutable Gaudi::Accumulators::StaticHistogram<1> m_residualHist{this, "track_hits_distance_closest_approach", "Track-hit Distances", {100, 0, 1, "Distance [mm];Entries"}};
 
 };
 

Tracking/components/TracksFromGenParticles.cpp

@@ -3,7 +3,7 @@
 // edm4hep
 #include "edm4hep/MCParticleCollection.h"
 #include "edm4hep/TrackCollection.h"
-#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/TrackMCParticleLinkCollection.h"
 
 // marlin
 #include <marlinutil/HelixClass_double.h>
@@ -27,7 +27,7 @@
  */
 
 struct TracksFromGenParticles final
-  : k4FWCore::MultiTransformer<std::tuple<edm4hep::TrackCollection, edm4hep::MCRecoTrackParticleAssociationCollection>(const edm4hep::MCParticleCollection&)> {
+  : k4FWCore::MultiTransformer<std::tuple<edm4hep::TrackCollection, edm4hep::TrackMCParticleLinkCollection>(const edm4hep::MCParticleCollection&)> {
   TracksFromGenParticles(const std::string& name, ISvcLocator* svcLoc)
       : MultiTransformer(
             name, svcLoc,
@@ -36,10 +36,10 @@ struct TracksFromGenParticles final
             KeyValues("OutputMCRecoTrackParticleAssociation", {"TracksFromGenParticlesAssociation"})}) {
   }
 
-std::tuple<edm4hep::TrackCollection, edm4hep::MCRecoTrackParticleAssociationCollection> operator()(const edm4hep::MCParticleCollection& genParticleColl) const override {
+std::tuple<edm4hep::TrackCollection, edm4hep::TrackMCParticleLinkCollection> operator()(const edm4hep::MCParticleCollection& genParticleColl) const override {
 
     auto outputTrackCollection = edm4hep::TrackCollection();
-    auto MCRecoTrackParticleAssociationCollection = edm4hep::MCRecoTrackParticleAssociationCollection();
+    auto MCRecoTrackParticleAssociationCollection = edm4hep::TrackMCParticleLinkCollection();
 
     for (const auto& genParticle : genParticleColl) {
       debug() << "Particle decayed in tracker: " << genParticle.isDecayedInTracker() << endmsg;
@@ -76,9 +76,9 @@ std::tuple<edm4hep::TrackCollection, edm4hep::MCRecoTrackParticleAssociationColl
       outputTrackCollection.push_back(trackFromGen);
 
       // Building the association between tracks and genParticles
-      auto MCRecoTrackParticleAssociation = edm4hep::MutableMCRecoTrackParticleAssociation();
-      MCRecoTrackParticleAssociation.setRec(trackFromGen);
-      MCRecoTrackParticleAssociation.setSim(genParticle);
+      auto MCRecoTrackParticleAssociation = edm4hep::MutableTrackMCParticleLink();
+      MCRecoTrackParticleAssociation.setFrom(trackFromGen);
+      MCRecoTrackParticleAssociation.setTo(genParticle);
       MCRecoTrackParticleAssociationCollection.push_back(MCRecoTrackParticleAssociation);
     }
     return std::make_tuple(std::move(outputTrackCollection), std::move(MCRecoTrackParticleAssociationCollection));

VTXdigi/include/VTXdigitizer.h

@@ -30,11 +30,13 @@ namespace edm4hep {
 
 #include "DDSegmentation/BitFieldCoder.h"
 
+#include <vector>
+
 /** @class VTXdigitizer
  *
  *  Algorithm for creating digitized (meaning 'reconstructed' for now) vertex detector hits (edm4hep::TrackerHit3D) from Geant4 hits (edm4hep::SimTrackerHit).
  *  
- *  @author Brieuc Francois
+ *  @author Brieuc Francois, Armin Ilg
  *  @date   2023-03
  *
  */
@@ -76,13 +78,13 @@ class VTXdigitizer : public Gaudi::Algorithm {
   dd4hep::VolumeManager m_volman;
 
   // x resolution in um
-  FloatProperty m_x_resolution{this, "xResolution", 0.1, "Spatial resolution in the x direction [um] (r-phi direction in barrel, z direction in disks)"};
+  Gaudi::Property<std::vector<float>> m_x_resolution{this, "xResolution", {0.1}, "Spatial resolutions in the x direction per layer [um] (r-phi direction in barrel, z direction in disks)"};
 
   // y resolution in um
-  FloatProperty m_y_resolution{this, "yResolution", 0.1, "Spatial resolution in the y direction [um] (r direction in barrel, r-phi direction in disks)"};
+  Gaudi::Property<std::vector<float>> m_y_resolution{this, "yResolution", {0.1}, "Spatial resolutions in the y direction per layer [um] (r direction in barrel, r-phi direction in disks)"};
 
   // t resolution in ns
-  FloatProperty m_t_resolution{this, "tResolution", 0.1, "Time resolution [ns]"};
+  Gaudi::Property<std::vector<float>> m_t_resolution{this, "tResolution", {0.1}, "Time resolutions per layer [ns]"};
 
   // Surface manager used to project hits onto sensitive surface with forceHitsOntoSurface argument
   mutable const dd4hep::rec::SurfaceMap* _map;
@@ -94,7 +96,7 @@ class VTXdigitizer : public Gaudi::Algorithm {
   IRndmGenSvc* m_randSvc;
 
   // Gaussian random number generator used for smearing
-  Rndm::Numbers m_gauss_x;
-  Rndm::Numbers m_gauss_y;
-  Rndm::Numbers m_gauss_time;
+  std::vector<Rndm::Numbers> m_gauss_x_vec;
+  std::vector<Rndm::Numbers> m_gauss_y_vec;
+  std::vector<Rndm::Numbers> m_gauss_t_vec;
 };

VTXdigi/src/VTXdigitizer.cpp

@@ -17,17 +17,29 @@ StatusCode VTXdigitizer::initialize() {
     error() << "Couldn't get RndmGenSvc!" << endmsg;
     return StatusCode::FAILURE;
   }
-  if (m_gauss_x.initialize(m_randSvc, Rndm::Gauss(0., m_x_resolution)).isFailure()) {
-    error() << "Couldn't initialize RndmGenSvc!" << endmsg;
-    return StatusCode::FAILURE;
+
+  m_gauss_x_vec.resize(m_x_resolution.size());
+  for (size_t i = 0; i < m_x_resolution.size(); ++i) {
+    if (m_gauss_x_vec[i].initialize(m_randSvc, Rndm::Gauss(0., m_x_resolution[i])).isFailure()) {
+      error() << "Couldn't initialize RndmGenSvc!" << endmsg;
+      return StatusCode::FAILURE;
+    }
   }
-  if (m_gauss_y.initialize(m_randSvc, Rndm::Gauss(0., m_y_resolution)).isFailure()) {
-    error() << "Couldn't initialize RndmGenSvc!" << endmsg;
-    return StatusCode::FAILURE;
+
+  m_gauss_y_vec.resize(m_y_resolution.size());
+  for (size_t i = 0; i < m_y_resolution.size(); ++i) {  
+    if (m_gauss_y_vec[i].initialize(m_randSvc, Rndm::Gauss(0., m_y_resolution[i])).isFailure()) {
+      error() << "Couldn't initialize RndmGenSvc!" << endmsg;
+      return StatusCode::FAILURE;
+    }
   }
-  if (m_gauss_time.initialize(m_randSvc, Rndm::Gauss(0., m_t_resolution)).isFailure()) {
-    error() << "Couldn't initialize RndmGenSvc!" << endmsg;
-    return StatusCode::FAILURE;
+
+  m_gauss_t_vec.resize(m_t_resolution.size());
+  for (size_t i = 0; i < m_t_resolution.size(); ++i) {  
+    if (m_gauss_t_vec[i].initialize(m_randSvc, Rndm::Gauss(0., m_t_resolution[i])).isFailure()) {
+      error() << "Couldn't initialize RndmGenSvc!" << endmsg;
+      return StatusCode::FAILURE;
+    }
   }
 
   // check if readout exists
@@ -39,6 +51,13 @@ StatusCode VTXdigitizer::initialize() {
   // set the cellID decoder
   m_decoder = m_geoSvc->getDetector()->readout(m_readoutName).idSpec().decoder(); // Can be used to access e.g. layer index: m_decoder->get(cellID, "layer"),
 
+  if (m_decoder->fieldDescription().find("layer") == std::string::npos){
+    error() 
+      << " Readout " << m_readoutName << " does not contain layer id!"
+      << endmsg;
+    return StatusCode::FAILURE;
+  }
+
   // retrieve the volume manager
   m_volman = m_geoSvc->getDetector()->volumeManager();
 
@@ -126,22 +145,21 @@ StatusCode VTXdigitizer::execute(const EventContext&) const {
 
     // Smear the hit in the local sensor coordinates
     double digiHitLocalPosition[3];
-    if (m_readoutName == "VTXIBCollection" ||
-        m_readoutName == "VTXOBCollection" ||
-        m_readoutName == "VertexBarrelCollection" ||
+    int iLayer = m_decoder->get(cellID, "layer");      
+    debug() << "readout: " << m_readoutName << ", layer id: " << iLayer << endmsg;
+    if (m_readoutName == "VertexBarrelCollection" ||
         m_readoutName == "SiWrBCollection") {  // In barrel, the sensor box is along y-z
       digiHitLocalPosition[0] = simHitLocalPositionVector.x();
-      digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_x.shoot() * dd4hep::mm;
-      digiHitLocalPosition[2] = simHitLocalPositionVector.z() + m_gauss_y.shoot() * dd4hep::mm;
-    } else if (m_readoutName == "VTXDCollection" ||
-               m_readoutName == "VertexEndcapCollection" ||
+      digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_x_vec[iLayer].shoot() * dd4hep::mm;
+      digiHitLocalPosition[2] = simHitLocalPositionVector.z() + m_gauss_y_vec[iLayer].shoot() * dd4hep::mm;
+    } else if (m_readoutName == "VertexEndcapCollection" ||
                m_readoutName == "SiWrDCollection") {  // In the disks, the sensor box is already in x-y
-      digiHitLocalPosition[0] = simHitLocalPositionVector.x() + m_gauss_x.shoot() * dd4hep::mm;
-      digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_y.shoot() * dd4hep::mm;
+      digiHitLocalPosition[0] = simHitLocalPositionVector.x() + m_gauss_x_vec[iLayer].shoot() * dd4hep::mm;
+      digiHitLocalPosition[1] = simHitLocalPositionVector.y() + m_gauss_y_vec[iLayer].shoot() * dd4hep::mm;
       digiHitLocalPosition[2] = simHitLocalPositionVector.z();
     } else {
       error()
-          << "VTX readout name (m_readoutName) unknown!"
+          << "VTX readout name (m_readoutName) unknown or xResolution/yResolution/tResolution not defining all detector layer resolutions!"
           << endmsg;
       return StatusCode::FAILURE;
     }
@@ -167,7 +185,7 @@ StatusCode VTXdigitizer::execute(const EventContext&) const {
     output_digi_hit.setPosition(digiHitGlobalPositionVector);
 
     // Apply time smearing
-    output_digi_hit.setTime(input_sim_hit.getTime() + m_gauss_time.shoot());
+    output_digi_hit.setTime(input_sim_hit.getTime() + m_gauss_t_vec[iLayer].shoot());
 
     output_digi_hit.setCellID(cellID);