Merge remote-tracking branch 'origin' into fix_CellID_FarClustering

.github/eicrecon.json

@@ -0,0 +1,15 @@
+{
+    "problemMatcher": [
+        {
+            "owner": "eicrecon",
+            "pattern": [
+                {
+                    "regexp": "^\\[(.*)\\]\\s\\[(error)\\]\\s*(.*)$",
+                    "file": 1,
+                    "severity": 2,
+                    "message": 3
+                }
+            ]
+        }
+    ]
+}

.github/ubsan.supp

@@ -1,3 +1,7 @@
 implicit-integer-sign-change:JsonMaterialDecorator
-implicit-integer-sign-change:/usr/local/include/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h
+implicit-integer-sign-change:/opt/local/include/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h
 implicit-integer-sign-change:TObject
+implicit-signed-integer-truncation:/opt/local/include/boost/iostreams/filter/gzip.hpp
+implicit-integer-sign-change:/opt/local/include/eigen3/Eigen/src/Core/arch/SSE/Complex.h
+implicit-integer-sign-change:/opt/local/include/root/TStorage.h
+implicit-unsigned-integer-truncation:/opt/local/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp

.pre-commit-config.yaml

@@ -7,7 +7,7 @@ repos:
   - id: end-of-file-fixer
   - id: trailing-whitespace
 - repo: https://github.com/codespell-project/codespell
-  rev: v2.2.6
+  rev: v2.3.0
   hooks:
   - id: codespell
 - repo: https://github.com/Lucas-C/pre-commit-hooks

CMakeLists.txt

@@ -31,7 +31,7 @@ cmake_minimum_required(VERSION 3.24)
 # 3.27: find_package() uses upper-case <PACKAGENAME>_ROOT variables
 cmake_policy(SET CMP0144 NEW)
 
-project(EICrecon)
+project(EICrecon LANGUAGES CXX)
 
 # CMake includes
 include(CheckCXXCompilerFlag)
@@ -197,8 +197,7 @@ endif()
 find_package(DD4hep REQUIRED)
 
 # ACTS cmake-lint: disable=C0103
-find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo
-                                      PluginDD4hep)
+find_package(Acts REQUIRED COMPONENTS Core PluginDD4hep PluginJson)
 set(Acts_VERSION_MIN "19.0.0")
 set(Acts_VERSION
     "${Acts_VERSION_MAJOR}.${Acts_VERSION_MINOR}.${Acts_VERSION_PATCH}")
@@ -233,6 +232,7 @@ endif()
 option(USE_ONNX "Compile with ONNX support" ON)
 if(${USE_ONNX})
   find_package(onnxruntime CONFIG)
+  add_compile_definitions(USE_ONNX)
 endif()
 
 # Add CMake additional functionality:

cmake/jana_plugin.cmake

@@ -269,8 +269,7 @@ endmacro()
 macro(plugin_add_acts _name)
 
   if(NOT Acts_FOUND)
-    find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo
-                                          PluginJson PluginDD4hep)
+    find_package(Acts REQUIRED COMPONENTS Core PluginDD4hep PluginJson)
     set(Acts_VERSION_MIN "20.2.0")
     set(Acts_VERSION
         "${Acts_VERSION_MAJOR}.${Acts_VERSION_MINOR}.${Acts_VERSION_PATCH}")
@@ -291,10 +290,8 @@ macro(plugin_add_acts _name)
   plugin_link_libraries(
     ${PLUGIN_NAME}
     ActsCore
-    ActsPluginIdentification
-    ActsPluginTGeo
-    ActsPluginJson
     ActsPluginDD4hep
+    ActsPluginJson
     ${ActsCore_PATH}/${CMAKE_SHARED_LIBRARY_PREFIX}ActsExamplesFramework${CMAKE_SHARED_LIBRARY_SUFFIX}
   )
   if(${_name}_WITH_LIBRARY)

docs/howto/visualize_callgraph.md

@@ -48,20 +48,20 @@ By default `eicrecon` activates the full reconstruction. This will
 result in a very busy callgraph that can be hard to read if you
 are only interested in one detector. (See [here](https://eic.github.io/EICrecon/howto/callgraphs/all.png).)
 To activate for a specific factory, you can specify the relevant output
-collection using the `podio:output_include_collections` parameter.
+collection using the `podio:output_collections` parameter.
 Here is an example for the *EcalEndcapN*:
 
 ~~~bash
 eicrecon -Pplugins=janadot \
-   -Ppodio:output_include_collections=EcalEndcapNMergedClusters \
+   -Ppodio:output_collections=EcalEndcapNMergedClusters \
    sim_file.edm4hep.root
 ~~~
 
 The above will produce something like the following:
 ![EcalEndcapNMergedClusters](callgraphs/EcalEndcapNMergedClusters.png)
 
 In the plot, the oval at the top is the JEventProcessor making the
-initial request for objects. (This is where the *podio:output_include_collections*
+initial request for objects. (This is where the *podio:output_collections*
 parameter is used). The rectangles in the middle represent factories (i.e. algorithms).
 The green trapezoid at the bottom is the object read from the input file.
 This is representative of the flow JANA implements. Namely, the requests for

docs/tutorial/02-work-environment.md

@@ -137,11 +137,11 @@ eicrecon 2022-09-26_ncdis10x100_minq2-1_100ev.edm4hep.root
 
 ## Generating a podio output file
 To write reconstructed values to an output file, you need to tell *eicrecon* what to write.
-There are several options available, but the mosrt useful one is *podio:output_include_collections*.
+There are several options available, but the mosrt useful one is *podio:output_collections*.
 This is a comma separated list of colelctions to write to the output file. For example:
 
 ```console
-eicrecon -Ppodio:output_include_collections=ReconstructedParticles 2022-09-26_ncdis10x100_minq2-1_100ev.edm4hep.root
+eicrecon -Ppodio:output_collections=ReconstructedParticles 2022-09-26_ncdis10x100_minq2-1_100ev.edm4hep.root
 ```
 
 To see a list of possible collections, run *eicrecon -L* .

docs/tutorial/04-factory.md

@@ -287,7 +287,7 @@ extern "C" {
 Finally, we go ahead and trigger the factory (remember, factories won't do anything unless activated by a JEventProcessor). You can open the
 
 ~~~ bash
-eicrecon in.root -Ppodio:output_file=out.root -Ppodio:output_include_collections=EcalEndcapNIslandClusters -Pjana:nevents=10
+eicrecon in.root -Ppodio:output_file=out.root -Ppodio:output_collections=EcalEndcapNIslandClusters -Pjana:nevents=10
 ~~~
 
 

docs/tutorial/05-contributing.md

@@ -52,7 +52,7 @@ More on the EIC contribution guide is in [this tutorial][eic-environment-tutoria
 
 ## Coding style
 
-One can find coding style and other contributins policies at [CONTRIBUTING.md](https://github.com/eic/EICrecon/blob/main/CONTRIBUTING.md). It is yet to be finished but one can find current decisions on coding style there
+One can find coding style and other contributions policies at [CONTRIBUTING.md](https://github.com/eic/EICrecon/blob/main/CONTRIBUTING.md). It is yet to be finished but one can find current decisions on coding style there
 
 
 ## References

src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc

@@ -161,14 +161,15 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
   // Used to optionally constrain the cluster eta to those of the contributing hits
   float minHitEta = std::numeric_limits<float>::max();
   float maxHitEta = std::numeric_limits<float>::min();
-  auto time       = pcl.getHits()[0].getTime();
-  auto timeError  = pcl.getHits()[0].getTimeError();
+  auto time      = 0;
+  auto timeError = 0;
   for (unsigned i = 0; i < pcl.getHits().size(); ++i) {
     const auto& hit   = pcl.getHits()[i];
     const auto weight = pcl.getWeights()[i];
     debug("hit energy = {} hit weight: {}", hit.getEnergy(), weight);
     auto energy = hit.getEnergy() * weight;
     totalE += energy;
+    time += (hit.getTime() - time) * energy / totalE;
     cl.addToHits(hit);
     cl.addToHitContributions(energy);
     const float eta = edm4hep::utils::eta(hit.getPosition());
@@ -228,12 +229,6 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
 
   // Additional convenience variables
 
-  // best estimate on the cluster direction is the cluster position
-  // for simple 2D CoG clustering
-  cl.setIntrinsicTheta(edm4hep::utils::anglePolar(cl.getPosition()));
-  cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(cl.getPosition()));
-  // TODO errors
-
   //_______________________________________
   // Calculate cluster profile:
   //    radius,
@@ -248,12 +243,11 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
   Eigen::Vector3f sum1_3D = Eigen::Vector3f::Zero();
   Eigen::Vector2cf eigenValues_2D = Eigen::Vector2cf::Zero();
   Eigen::Vector3cf eigenValues_3D = Eigen::Vector3cf::Zero();
-  //the axis is the direction of the eigenvalue corresponding to the largest eigenvalue.
-  double axis_x=0, axis_y=0, axis_z=0;
-  if (cl.getNhits() > 1) {
+  // the axis is the direction of the eigenvalue corresponding to the largest eigenvalue.
+  edm4hep::Vector3f axis;
 
+  if (cl.getNhits() > 1) {
     for (const auto& hit : pcl.getHits()) {
-
       float w = weightFunc(hit.getEnergy(), totalE, logWeightBase, 0);
 
       // theta, phi
@@ -306,19 +300,15 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
             return std::real(a) < std::real(b);
         }
       );
-      auto axis = eigenvectors.col(std::distance(
+      auto axis_eigen = eigenvectors.col(std::distance(
             eigenValues_3D.begin(),
             max_eigenvalue_it
         ));
-      axis_x=axis(0,0).real();
-      axis_y=axis(1,0).real();
-      axis_z=axis(2,0).real();
-      double norm=sqrt(axis_x*axis_x+axis_y*axis_y+axis_z*axis_z);
-      if (norm!=0){
-        axis_x/=norm;
-        axis_y/=norm;
-        axis_z/=norm;
-      }
+      axis = {
+        axis_eigen(0,0).real(),
+        axis_eigen(1,0).real(),
+        axis_eigen(2,0).real(),
+      };
     }
   }
 
@@ -329,10 +319,22 @@ std::optional<edm4eic::MutableCluster> CalorimeterClusterRecoCoG::reconstruct(co
   cl.addToShapeParameters( eigenValues_3D[0].real() ); // 3D x-y-z cluster width 1
   cl.addToShapeParameters( eigenValues_3D[1].real() ); // 3D x-y-z cluster width 2
   cl.addToShapeParameters( eigenValues_3D[2].real() ); // 3D x-y-z cluster width 3
-  //last 3 shape parameters are the components of the axis direction
-  cl.addToShapeParameters( axis_x );
-  cl.addToShapeParameters( axis_y );
-  cl.addToShapeParameters( axis_z );
+
+
+  double dot_product = cl.getPosition() * axis;
+  if (dot_product < 0) {
+    axis = -1 * axis;
+  }
+
+  if (m_cfg.longitudinalShowerInfoAvailable) {
+    cl.setIntrinsicTheta(edm4hep::utils::anglePolar(axis));
+    cl.setIntrinsicPhi(edm4hep::utils::angleAzimuthal(axis));
+    // TODO intrinsicDirectionError
+  } else {
+    cl.setIntrinsicTheta(NAN);
+    cl.setIntrinsicPhi(NAN);
+  }
+
   return std::move(cl);
 }
 

src/algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h

@@ -23,6 +23,8 @@ namespace eicrecon {
         std::vector<double> logWeightBaseCoeffs{};
         double logWeightBase_Eref = 50 * dd4hep::MeV;
 
+        bool longitudinalShowerInfoAvailable = false;
+
         // Constrain the cluster position eta to be within
         // the eta of the contributing hits. This is useful to avoid edge effects
         // for endcaps.

src/algorithms/calorimetry/CalorimeterHitDigi.cc

@@ -23,6 +23,7 @@
 #include <edm4hep/CaloHitContributionCollection.h>
 #include <fmt/core.h>
 #include <podio/RelationRange.h>
+#include <algorithm>
 #include <cmath>
 #include <cstddef>
 #include <gsl/pointers>
@@ -181,10 +182,13 @@ void CalorimeterHitDigi::process(
                      std::pow(m_cfg.eRes[2] / (edep), 2)
                   )
                 : 0;
-        double    corrMeanScale_value = corrMeanScale(leading_hit);
-        double    ped     = m_cfg.pedMeanADC + m_gaussian(m_generator) * m_cfg.pedSigmaADC;
-        unsigned long long adc     = std::llround(ped + edep * corrMeanScale_value * ( 1.0 + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC);
-        unsigned long long tdc     = std::llround((time + m_gaussian(m_generator) * tRes) * stepTDC);
+        double corrMeanScale_value = corrMeanScale(leading_hit);
+
+        double ped = m_cfg.pedMeanADC + m_gaussian(m_generator) * m_cfg.pedSigmaADC;
+
+        // Note: both adc and tdc values must be positive numbers to avoid integer wraparound
+        unsigned long long adc = std::max(std::llround(ped + edep * corrMeanScale_value * (1.0 + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC), 0LL);
+        unsigned long long tdc = std::llround((time + m_gaussian(m_generator) * tRes) * stepTDC);
 
         if (edep> 1.e-3) trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t corrMeanScale: {}", edep, adc, time, m_cfg.capTime, tdc, corrMeanScale_value);
         rawhits->create(

src/algorithms/calorimetry/CalorimeterHitReco.cc

@@ -21,6 +21,7 @@
 #include <Math/GenVector/Cartesian3D.h>
 #include <Math/GenVector/DisplacementVector3D.h>
 #include <algorithms/service.h>
+#include <edm4eic/EDM4eicVersion.h>
 #include <fmt/core.h>
 #include <fmt/format.h>
 #include <algorithm>
@@ -172,6 +173,11 @@ void CalorimeterHitReco::process(
             continue;
         }
 
+        if (rh.getAmplitude() > m_cfg.capADC) {
+            error("Encountered hit with amplitude {} outside of ADC capacity {}", rh.getAmplitude(), m_cfg.capADC);
+            continue;
+        }
+
         // get layer and sector ID
         const int lid =
                 id_dec != nullptr && !m_cfg.layerField.empty() ? static_cast<int>(id_dec->get(cellID, layer_idx)) : -1;
@@ -265,6 +271,9 @@ void CalorimeterHitReco::process(
         const decltype(edm4eic::CalorimeterHitData::local) local_position(pos.x() / dd4hep::mm, pos.y() / dd4hep::mm,
                                                                        pos.z() / dd4hep::mm);
 
+#if EDM4EIC_VERSION_MAJOR >= 7
+        auto recohit =
+#endif
         recohits->create(
             rh.getCellID(),
             energy,
@@ -276,6 +285,9 @@ void CalorimeterHitReco::process(
             sid,
             lid,
             local_position);
+#if EDM4EIC_VERSION_MAJOR >= 7
+        recohit.setRawHit(rh);
+#endif
     }
 }