From 1793804701fa43114f25ce0e1984ac62b6d201d4 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Thu, 28 Sep 2023 15:36:38 +0200
Subject: [PATCH 01/24] Add a DDPlanarDigiProcessor

---
 k4Reco/CMakeLists.txt                         |  13 +
 .../components/DDPlanarDigiProcessor.cpp      | 358 ++++++++++++++++++
 .../components/DDPlanarDigiProcessor.h        | 114 ++++++
 .../options/runDDPlanarDigiProcessor.py       |  59 +++
 4 files changed, 544 insertions(+)
 create mode 100644 k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
 create mode 100644 k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
 create mode 100644 k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py

diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt
index 04c18f6..618b309 100644
--- a/k4Reco/CMakeLists.txt
+++ b/k4Reco/CMakeLists.txt
@@ -16,3 +16,16 @@ WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 ]]
+
+set(sources DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+)
+gaudi_add_module(k4RecoDDPlanarDigi
+                 SOURCES ${sources}
+                 LINK Gaudi::GaudiKernel
+                      Gaudi::GaudiAlgLib
+                      k4FWCore::k4FWCore
+                      k4FWCore::k4Interface
+                      EDM4HEP::edm4hep
+                      DD4hep::DDRec
+                      DD4hep::DDCore
+                      )
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
new file mode 100644
index 0000000..70bc200
--- /dev/null
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
@@ -0,0 +1,358 @@
+/*
+ * Copyright (c) 2014-2023 Key4hep-Project.
+ *
+ * This file is part of Key4hep.
+ * See https://key4hep.github.io/key4hep-doc/ for further info.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#include "DDPlanarDigiProcessor.h"
+
+#include "edm4hep/SimTrackerHit.h"
+#include "edm4hep/TrackerHitPlaneCollection.h"
+#include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
+
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/BitFieldCoder.h"
+
+#include <TMath.h>
+#include <TFile.h>
+
+#include <cmath>
+#include <sstream>
+#include <iostream>
+
+enum {
+  hu = 0,
+  hv,
+  hT,
+  hitE,
+  hitsAccepted,
+  diffu,
+  diffv,
+  diffT,
+  hSize 
+};
+
+DDPlanarDigiProcessor::DDPlanarDigiProcessor(const std::string& name, ISvcLocator* svcLoc)
+    : MultiTransformer(
+          name, svcLoc,
+          {
+           KeyValue("SimTrackerHitCollectionName", "SimTrackerHitss"),
+           },
+          {KeyValue("TrackerHitCollectionName", "VTXTrackerHits"),
+           KeyValue("SimTrkHitRelCollection", "VTXTrackerHitRelations")}) {
+
+
+  auto randSvc = service<IRndmGenSvc>("RndmGenSvc", true);
+  if (!randSvc || !rng.initialize(randSvc, Rndm::Gauss(0, 1))) {
+    error() << "Unable to create Random generator" << endmsg;
+  }
+
+  m_uidSvc = service<IUniqueIDGenSvc>("UniqueIDGenSvc", true);
+  if (!m_uidSvc) {
+    error() << "Unable to get UniqueIDGenSvc" << endmsg;
+  }
+
+  if(m_resULayer.size() != m_resVLayer.size()) {
+    error() << "DDPlanarDigiProcessor - Inconsistent number of resolutions given for U and V coordinate: " 
+       << "ResolutionU  :" <<   m_resULayer.size() << " != ResolutionV : " <<  m_resVLayer.size();
+
+    throw std::runtime_error("DDPlanarDigiProcessor: Inconsistent number of resolutions given for U and V coordinate");
+  }
+
+  m_histograms.resize(hSize);
+
+  m_histograms[hu] = new TH1F("hu", "smearing u", 50, -5., +5.);
+  m_histograms[hv] = new TH1F("hv", "smearing v", 50, -5., +5.);
+  m_histograms[hT] = new TH1F("hT", "smearing time", 50, -5., +5.);
+
+  m_histograms[diffu] = new TH1F("diffu", "diff u", 1000, -5., +5.);
+  m_histograms[diffv] = new TH1F("diffv", "diff v", 1000, -5., +5.);
+  m_histograms[diffT] = new TH1F("diffT", "diff time", 1000, -5., +5.);
+
+  m_histograms[hitE] = new TH1F("hitE", "hitEnergy in keV", 1000, 0, 200);
+  m_histograms[hitsAccepted] = new TH1F("hitsAccepted", "Fraction of accepted hits [%]", 201, 0, 100.5);
+
+  m_geoSvc = serviceLocator()->service(m_geoSvcName);
+}
+
+StatusCode
+DDPlanarDigiProcessor::initialize() {
+  const auto detector = m_geoSvc->getDetector();
+
+  const auto surfMan = detector->extension<dd4hep::rec::SurfaceManager>();
+  dd4hep::DetElement det = detector->detector(m_subDetName.value());
+  surfaceMap = surfMan->map(m_subDetName.value());
+
+  if(!surfaceMap) {   
+    std::stringstream err ; err << " Could not find surface map for detector: " 
+                                << det.name() << " in SurfaceManager ";
+    throw std::runtime_error(err.str());
+  }
+  return StatusCode::SUCCESS;
+}
+
+
+std::tuple<TrackerHitPlaneColl, Association>
+DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits) const {
+
+  // TODO: Set seed
+  // gslrng_set( rng, Global::EVENTSEEDER->getSeed(this) );   
+  // debug() << "seed set to " << Global::EVENTSEEDER->getSeed(this) << endmsg;
+  auto seed = m_uidSvc->getUniqueID(1, 2, "hello");
+  
+  int nCreatedHits=0;
+  int nDismissedHits=0;
+
+  auto trkhitVec = edm4hep::TrackerHitPlaneCollection();
+  auto thsthcol = edm4hep::MCRecoTrackerHitPlaneAssociationCollection();
+
+  std::string cellIDEncodingString = m_geoSvc->constantAsString(m_encodingStringVariable.value());
+  dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(cellIDEncodingString); 
+  
+  int nSimHits = simTrackerHits.size();
+  std::string collName;
+  (void)this->getProperty("SimTrackerHitCollectionName", collName);
+  debug() << "Processing collection " << collName << " with " << simTrackerHits.size() << " hits ... " << endmsg;
+
+  for(const auto& hit : simTrackerHits) {
+    m_histograms[hitE]->Fill( hit.getEDep() * (dd4hep::GeV / dd4hep::keV) );
+
+    if(hit.getEDep() < m_minEnergy) {
+      debug() << "Hit with insufficient energy " << hit.getEDep() * (dd4hep::GeV / dd4hep::keV) << " keV" << endmsg;
+      continue;
+    }
+
+    const int cellID0 = hit.getCellID();
+
+    // get the measurement surface for this hit using the CellID
+    dd4hep::rec::SurfaceMap::const_iterator sI = surfaceMap->find(cellID0);
+
+    if(sI == surfaceMap->end()){
+      std::stringstream err; err << " DDPlanarDigiProcessor::processEvent(): no surface found for cellID : " << std::endl;
+                                  // <<   cellid_decoder( hit ).valueString() ;
+      throw std::runtime_error(err.str());
+    }
+
+    const dd4hep::rec::ISurface* surf = sI->second;
+    int layer = bitFieldCoder.get(cellID0, "layer");
+
+    dd4hep::rec::Vector3D oldPos(hit.getPosition()[0], hit.getPosition()[1], hit.getPosition()[2]);
+    dd4hep::rec::Vector3D newPos;
+
+    //  Check if Hit is inside sensitive 
+    if (!surf->insideBounds(dd4hep::mm * oldPos)) {
+      
+      // debug() << "  hit at " << oldPos 
+      //         << " " << cellid_decoder( hit).valueString() 
+      //         << " is not on surface " 
+      //         << *surf  
+      //         << " distance: " << surf->distance(  dd4hep::mm * oldPos )
+      //         << endmsg;        
+
+      if(m_forceHitsOntoSurface){
+        dd4hep::rec::Vector2D lv = surf->globalToLocal(dd4hep::mm * oldPos);
+        dd4hep::rec::Vector3D oldPosOnSurf = (1./dd4hep::mm) * surf->localToGlobal(lv); 
+        
+        debug() << " moved to " << oldPosOnSurf << " distance " << (oldPosOnSurf-oldPos).r()
+                << endmsg;        
+          
+        oldPos = oldPosOnSurf;
+
+      } else {
+        ++nDismissedHits;
+        continue; 
+      }
+    }
+
+    // Smear time of the hit and apply the time window cut if needed
+    double hitT = hit.getTime();
+    
+    if (m_resTLayer.size() and m_resTLayer[0] > 0) {
+      float resT = m_resTLayer.size() > 1 ? m_resTLayer[layer] : m_resTLayer[0];
+
+      // For this and other rng calls, use the fact that drawing from a
+      // Gaussing with mean mu and variance sigma^2 is the same as drawing from
+      // a normal distribution and then multiplying by sigma and adding mu
+      double tSmear  = resT > 0 ? rng() * resT : 0;
+      m_histograms[hT]->Fill(resT > 0 ? tSmear / resT : 0);
+      m_histograms[diffT]->Fill(tSmear);
+
+      hitT += tSmear;
+      debug() << "smeared hit at T: " << hit.getTime() << " ns to T: " << hitT << " ns according to resolution: " << resT << " ns" << endmsg;
+    }
+    
+    // Correcting for the propagation time
+    if (m_correctTimesForPropagation) {
+      double dt = oldPos.r() / (TMath::C()/1e6);
+      hitT -= dt;
+      debug() << "corrected hit at R: " << oldPos.r() << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << endmsg;
+    }
+    
+    // Skip the hit if its time is outside the acceptance time window
+    if (m_useTimeWindow) {
+      // TODO: Check that the length of the time window is OK
+      float timeWindow_min = m_timeWindowMin.size() > 1 ? m_timeWindowMin[layer] : m_timeWindowMin[0];
+      float timeWindow_max = m_timeWindowMax.size() > 1 ? m_timeWindowMax[layer] : m_timeWindowMax[0];
+      if ( hitT < timeWindow_min || hitT > timeWindow_max ) {
+        debug() << "hit at T: " << hit.getTime() << " smeared to: " << hitT << " is outside the time window: hit dropped"  << endmsg;
+        ++nDismissedHits;
+        continue;
+      }
+    }
+    
+    // Try to smear the hit position but ensure the hit is inside the sensitive region
+    dd4hep::rec::Vector3D u = surf->u();
+    dd4hep::rec::Vector3D v = surf->v();
+
+    // get local coordinates on surface
+    dd4hep::rec::Vector2D lv = surf->globalToLocal(dd4hep::mm * oldPos);
+    double uL = lv[0] / dd4hep::mm;
+    double vL = lv[1] / dd4hep::mm;
+
+    bool acceptHit = false;
+    int tries =  0;              
+    
+    // TODO: check lengths
+    float resU = m_resULayer.size() > 1 ? m_resULayer[layer] : m_resULayer[0];
+    float resV = m_resVLayer.size() > 1 ? m_resVLayer[layer] : m_resVLayer[0]; 
+
+    while(tries < m_maxTries) {
+      
+      // if( tries > 0 ) debug() << "retry smearing for " <<  cellid_decoder( hit ).valueString() << " : retries " << tries << endmsg;
+      
+      double uSmear  = rng() * resU;
+      double vSmear  = rng() * resV;
+      
+      dd4hep::rec::Vector3D newPosTmp;
+      if (m_isStrip){
+          if (m_subDetName == "SET"){
+            double xStripPos, yStripPos, zStripPos;
+            //Find intersection of the strip with the z=centerOfSensor plane to set it as the center of the SET strip
+            dd4hep::rec::Vector3D simHitPosSmeared = (1./dd4hep::mm) * ( surf->localToGlobal( dd4hep::rec::Vector2D( (uL+uSmear)*dd4hep::mm, 0.) ) );
+            zStripPos = surf->origin()[2] / dd4hep::mm;
+            double lineParam = (zStripPos - simHitPosSmeared[2])/v[2];
+            xStripPos = simHitPosSmeared[0] + lineParam*v[0];
+            yStripPos = simHitPosSmeared[1] + lineParam*v[1];
+            newPosTmp = dd4hep::rec::Vector3D(xStripPos, yStripPos, zStripPos);    
+          }
+          else{
+            newPosTmp = (1./dd4hep::mm) * ( surf->localToGlobal( dd4hep::rec::Vector2D( (uL+uSmear)*dd4hep::mm, 0. ) ) );    
+          }
+      }
+      else{
+        newPosTmp = (1./dd4hep::mm) * ( surf->localToGlobal( dd4hep::rec::Vector2D( (uL+uSmear)*dd4hep::mm, (vL+vSmear)*dd4hep::mm ) ) );
+      }
+
+      debug() << " hit at    : " << oldPos 
+              << " smeared to: " << newPosTmp
+              << " uL: " << uL 
+              << " vL: " << vL 
+              << " uSmear: " << uSmear
+              << " vSmear: " << vSmear
+              << endmsg;
+
+
+      if ( surf->insideBounds( dd4hep::mm * newPosTmp ) ) {    
+        acceptHit = true;
+        newPos     = newPosTmp;
+
+        m_histograms[hu]->Fill(  uSmear / resU ); 
+        m_histograms[hv]->Fill(  vSmear / resV ); 
+
+        m_histograms[diffu]->Fill( uSmear );
+        m_histograms[diffv]->Fill( vSmear );
+
+        break;  
+      }
+
+      // debug() << "  hit at " << newPosTmp 
+      //         << " " << cellid_decoder( hit).valueString() 
+      //         << " is not on surface " 
+      //         << " distance: " << surf->distance( dd4hep::mm * newPosTmp ) 
+      //         << endmsg;        
+      
+      ++tries;
+    }
+
+   if(!acceptHit) {
+     debug() << "hit could not be smeared within ladder after " << m_maxTries << "  tries: hit dropped"  << endmsg;
+     ++nDismissedHits;
+     continue; 
+   } 
+   
+   auto trkHit = trkhitVec.create();
+               
+   trkHit.setCellID( cellID0 );
+   
+   trkHit.setPosition( newPos.const_array()  );
+   trkHit.setTime( hitT );
+   trkHit.setEDep( hit.getEDep() );
+
+   float u_direction[2];
+   u_direction[0] = u.theta();
+   u_direction[1] = u.phi();
+   
+   float v_direction[2];
+   v_direction[0] = v.theta();
+   v_direction[1] = v.phi();
+
+   debug()  << " U[0] = "<< u_direction[0] << " U[1] = "<< u_direction[1] 
+                          << " V[0] = "<< v_direction[0] << " V[1] = "<< v_direction[1]
+                          << endmsg;
+
+   trkHit.setU(u_direction);
+   trkHit.setV(v_direction);
+   trkHit.setDu(resU);
+
+   if(m_isStrip) {
+     // store the resolution from the length of the wafer - in case a fitter might want to treat this as 2d hit ....
+     double stripRes = surf->length_along_v() / dd4hep::mm / std::sqrt(12);
+     trkHit.setDv(stripRes); 
+     // TODO: Set type?
+     // trkHit.setType( UTIL::set_bit( trkHit.getType() ,  UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ) );
+
+   } else {
+     trkHit.setDv(resV);
+   }
+
+   auto association = thsthcol.create();
+   association.setSim(hit);
+   association.setRec(trkHit);
+   
+   ++nCreatedHits;
+   
+  }
+
+  // Filling the fraction of accepted hits in the event
+  float accFraction = nSimHits > 0 ? float(nCreatedHits) / float(nSimHits) * 100.0 : 0.0;
+  m_histograms[hitsAccepted]->Fill(accFraction);
+  
+  debug() << "Created " << nCreatedHits << " hits, " << nDismissedHits << " hits  dismissed" << endmsg;
+  
+  return std::make_tuple(std::move(trkhitVec), std::move(thsthcol));
+}
+
+
+StatusCode
+DDPlanarDigiProcessor::finalize() {
+  auto file = TFile::Open("DDPlanarDigiProcessor.root", "RECREATE");
+  for (auto& h : m_histograms) {
+    h->Write();
+    delete h;
+  }
+  file->Close();
+  return StatusCode::SUCCESS;
+}
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
new file mode 100644
index 0000000..d907bb5
--- /dev/null
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
@@ -0,0 +1,114 @@
+/*
+ * Copyright (c) 2014-2023 Key4hep-Project.
+ *
+ * This file is part of Key4hep.
+ * See https://key4hep.github.io/key4hep-doc/ for further info.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+#ifndef DDPLANARDIGIPROCESSOR_H
+#define DDPLANARDIGIPROCESSOR_H
+
+#include "Gaudi/Property.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "k4Interface/IGeoSvc.h"
+#include "k4FWCore/BaseClass.h"
+#include "k4Interface/IUniqueIDGenSvc.h"
+
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitPlaneCollection.h"
+#include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
+
+#include "DDRec/SurfaceManager.h"
+
+#include <string>
+#include <TH1F.h>
+
+/** ======= DDPlanarDigiProcessor ========== <br>
+ * Creates TrackerHits from SimTrackerHits, smearing them according to the input parameters. 
+ * The positions of "digitized" TrackerHits are obtained by gaussian smearing positions
+ * of SimTrackerHits perpendicular and along the ladder according to the specified point resolutions. 
+ * The geometry of the surface is retreived from DDRec::Surface associated to the hit via cellID.
+ * 
+ * 
+ * <h4>Input collections and prerequisites</h4> 
+ * Processor requires a collection of SimTrackerHits <br>
+ * <h4>Output</h4>
+ * Processor produces collection of smeared TrackerHits<br>
+ * @param SimTrackHitCollectionName The name of input collection of SimTrackerHits <br>
+ * (default name VXDCollection) <br>
+ * @param TrackerHitCollectionName The name of output collection of smeared TrackerHits <br>
+ * (default name VTXTrackerHits) <br>
+ * @param ResolutionU resolution in direction of u (in mm) <br>
+ * (default value 0.004) <br>
+ * @param ResolutionV Resolution in direction of v (in mm) <br>
+ * (default value 0.004) <br>
+ * @param IsStrip whether the hits are 1 dimensional strip measurements <br>
+ * (default value false)
+ * @param Ladder_Number_encoded_in_cellID ladder number has been encoded in the cellID <br>
+ * (default value false) <br>
+ * @param Sub_Detector_ID ID of Sub-Detector using UTIL/ILDConf.h from lcio <br>
+ * (default value lcio::ILDDetID::VXD) <br>
+ * <br>
+ * 
+ * @author F.Gaede CERN/DESY, S. Aplin DESY
+ * @date Dec 2014
+ * Originally in https://github.com/iLCSoft/MarlinTrkProcessors/blob/master/source/Digitisers/include/DDPlanarDigiProcessor.h
+ */
+
+using SimTrackerHitCollection = edm4hep::SimTrackerHitCollection;
+
+using TrackerHitPlaneColl = edm4hep::TrackerHitPlaneCollection;
+using Association = edm4hep::MCRecoTrackerHitPlaneAssociationCollection;
+
+struct DDPlanarDigiProcessor final
+    : Gaudi::Functional::MultiTransformer<std::tuple<TrackerHitPlaneColl, Association>(const SimTrackerHitCollection&), BaseClass_t> {
+
+  DDPlanarDigiProcessor(const std::string& name, ISvcLocator* svcLoc);
+  
+  StatusCode initialize() override;
+  StatusCode finalize() override;
+
+  std::tuple<TrackerHitPlaneColl, Association> operator()(const SimTrackerHitCollection& simTrackerHits) const override;
+
+ private:
+  Gaudi::Property<std::string> m_subDetName{this, "SubDetectorName", "VXD", "Name of the subdetector"};
+  Gaudi::Property<bool> m_isStrip{this, "IsStrip", false, "Whether the hits are 1D strip hits"};
+  Gaudi::Property<std::vector<float>> m_resULayer{this, "ResolutionU", {0.004}, "Resolution in the direction of u; either one per layer or one for all layers"};
+  Gaudi::Property<std::vector<float>> m_resVLayer{this, "ResolutionV", {0.004}, "Resolution in the direction of v; either one per layer or one for all layers"};
+  Gaudi::Property<std::vector<float>> m_resTLayer{this, "ResolutionT", {0.004}, "Resolution in the direction of t; either one per layer or one for all layers. If the single entry is negative, disable time smearing. "};
+  Gaudi::Property<bool> m_forceHitsOntoSurface{this, "ForceHitsOntoSurface", false, "Project hits onto the surfoce in case they are not yet on the surface"};
+  Gaudi::Property<double> m_minEnergy{this, "MinEnergy", 0.0, "Minimum energy (GeV) of SimTrackerHit to be digitized"};
+
+  Gaudi::Property<bool> m_useTimeWindow{this, "UseTimeWindow", false, "Only accept hits with time (after smearing) within the specified time window (default: false)"};
+  Gaudi::Property<bool> m_correctTimesForPropagation{this, "CorrectTimesForPropagation", false, "Correct hit time for the propagation: radial distance/c (default: false)"};
+  Gaudi::Property<std::vector<float>> m_timeWindowMin{this, "TimeWindowMin", {-1e9}, "Minimum time (ns) of SimTrackerHit to be digitized"};
+  Gaudi::Property<std::vector<float>> m_timeWindowMax{this, "TimeWindowMax", {1e9}, "Maximum time (ns) of SimTrackerHit to be digitized"};
+  Gaudi::Property<std::string> m_encodingStringVariable{
+      this, "EncodingStringParameterName", "GlobalTrackerReadoutID",
+      "The name of the DD4hep constant that contains the Encoding string for tracking detectors"};
+  Gaudi::Property<std::string> m_geoSvcName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"};
+  Gaudi::Property<int> m_maxTries{this, "MaxTries", 10, "Maximum number of tries to find a valid surface for a hit"};
+
+  Rndm::Numbers rng;
+  const dd4hep::rec::SurfaceMap* surfaceMap;
+  std::vector<TH1F*> m_histograms;
+  SmartIF<IGeoSvc> m_geoSvc;
+  SmartIF<IUniqueIDGenSvc> m_uidSvc;
+  
+};
+
+DECLARE_COMPONENT(DDPlanarDigiProcessor)
+
+#endif
diff --git a/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py b/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
new file mode 100644
index 0000000..8ac0250
--- /dev/null
+++ b/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
@@ -0,0 +1,59 @@
+#
+# Copyright (c) 2020-2023 Key4hep-Project.
+#
+# This file is part of Key4hep.
+# See https://key4hep.github.io/key4hep-doc/ for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+from Gaudi.Configuration import INFO
+from Configurables import ApplicationMgr, k4DataSvc, PodioInput, PodioOutput
+from Configurables import DDPlanarDigiProcessor
+from Configurables import GeoSvc
+from Configurables import UniqueIDGenSvc
+import os
+
+id_service = UniqueIDGenSvc("UniqueIDGenSvc")
+
+geoservice = GeoSvc("GeoSvc")
+geoservice.detectors = [os.environ["K4GEO"]+"/CLIC/compact/CLIC_o3_v15/CLIC_o3_v15.xml"]
+geoservice.OutputLevel = INFO
+geoservice.EnableGeant4Geo = False
+
+processor = DDPlanarDigiProcessor()
+processor.SubDetectorName = "Vertex"
+processor.IsStrip = False
+processor.ResolutionU = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
+processor.ResolutionV = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
+processor.SimTrackerHitCollectionName = "VertexBarrelCollection"
+processor.SimTrkHitRelCollection = "VXDTrackerHitRelations"
+processor.TrackerHitCollectionName = "VXDTrackerHits"
+
+data_svc = k4DataSvc("EventDataSvc")
+data_svc.input = "input.root"
+
+inp = PodioInput()
+inp.collections = [
+    "VertexBarrelCollection",
+]
+
+out = PodioOutput("out")
+out.filename = "ttbar_edm4hep_digi.root"
+
+ApplicationMgr(TopAlg=[inp, processor, out],
+               EvtSel="NONE",
+               EvtMax=-1,
+               ExtSvc=[data_svc],
+               OutputLevel=INFO,
+               )
+

From df3362868b00bc2190916a1aff5ccfa2593cd978 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Thu, 28 Sep 2023 16:17:43 +0200
Subject: [PATCH 02/24] Add k4RecoPlugins to the CMakeLists

---
 k4Reco/CMakeLists.txt | 26 +++++++++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt
index 618b309..489b9ec 100644
--- a/k4Reco/CMakeLists.txt
+++ b/k4Reco/CMakeLists.txt
@@ -18,8 +18,9 @@ limitations under the License.
 ]]
 
 set(sources DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+            components/InitializeDD4hep.cpp
 )
-gaudi_add_module(k4RecoDDPlanarDigi
+gaudi_add_module(k4RecoPlugins
                  SOURCES ${sources}
                  LINK Gaudi::GaudiKernel
                       Gaudi::GaudiAlgLib
@@ -29,3 +30,26 @@ gaudi_add_module(k4RecoDDPlanarDigi
                       DD4hep::DDRec
                       DD4hep::DDCore
                       )
+target_include_directories(k4RecoPlugins PUBLIC
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/include>
+  $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+install(TARGETS k4RecoPlugins
+  EXPORT ${CMAKE_PROJECT_NAME}Targets
+  RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+  LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+  PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/@{CMAKE_PROJECT_NAME}"
+  COMPONENT dev)
+
+include(CTest)
+
+set(GAUDI_GENCONF_DIR "genConfDir")
+
+function(set_test_env _testname)
+  set_property(TEST ${_testname} APPEND PROPERTY ENVIRONMENT
+    LD_LIBRARY_PATH=${CMAKE_BINARY_DIR}:$<TARGET_FILE_DIR:k4RecoPlugins>:$<TARGET_FILE_DIR:ROOT::Core>:$<TARGET_FILE_DIR:k4FWCore::k4FWCore>:$<TARGET_FILE_DIR:EDM4HEP::edm4hep>:$<TARGET_FILE_DIR:podio::podio>:$ENV{LD_LIBRARY_PATH}
+    PYTHONPATH=${CMAKE_BINARY_DIR}/${CMAKE_PROJECT_NAME}/${GAUDI_GENCONF_DIR}:$<TARGET_FILE_DIR:k4FWCore::k4FWCore>/../python:$ENV{PYTHONPATH}
+    PATH=$<TARGET_FILE_DIR:k4FWCore::k4FWCore>/../bin:$ENV{PATH}
+    K4PROJECTTEMPLATE=${CMAKE_CURRENT_LIST_DIR}/
+    )
+endfunction()

From 6aaf55d78b5579c6eb784c9207e142b7d49e8b82 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Thu, 28 Sep 2023 16:18:30 +0200
Subject: [PATCH 03/24] Remove an outdated line

---
 k4Reco/CMakeLists.txt | 1 -
 1 file changed, 1 deletion(-)

diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt
index 489b9ec..35fd9c9 100644
--- a/k4Reco/CMakeLists.txt
+++ b/k4Reco/CMakeLists.txt
@@ -18,7 +18,6 @@ limitations under the License.
 ]]
 
 set(sources DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
-            components/InitializeDD4hep.cpp
 )
 gaudi_add_module(k4RecoPlugins
                  SOURCES ${sources}

From abfbc16e7c0470875c12c6f88287303c07de4471 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Tue, 3 Oct 2023 09:56:01 +0200
Subject: [PATCH 04/24] Use a random number generator from <random>

---
 .../components/DDPlanarDigiProcessor.cpp           | 14 ++++++--------
 .../components/DDPlanarDigiProcessor.h             |  3 ++-
 2 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
index 70bc200..c66b622 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
@@ -30,6 +30,7 @@
 #include <TFile.h>
 
 #include <cmath>
+#include <random>
 #include <sstream>
 #include <iostream>
 
@@ -55,11 +56,6 @@ DDPlanarDigiProcessor::DDPlanarDigiProcessor(const std::string& name, ISvcLocato
            KeyValue("SimTrkHitRelCollection", "VTXTrackerHitRelations")}) {
 
 
-  auto randSvc = service<IRndmGenSvc>("RndmGenSvc", true);
-  if (!randSvc || !rng.initialize(randSvc, Rndm::Gauss(0, 1))) {
-    error() << "Unable to create Random generator" << endmsg;
-  }
-
   m_uidSvc = service<IUniqueIDGenSvc>("UniqueIDGenSvc", true);
   if (!m_uidSvc) {
     error() << "Unable to get UniqueIDGenSvc" << endmsg;
@@ -112,6 +108,8 @@ DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits)
   // gslrng_set( rng, Global::EVENTSEEDER->getSeed(this) );   
   // debug() << "seed set to " << Global::EVENTSEEDER->getSeed(this) << endmsg;
   auto seed = m_uidSvc->getUniqueID(1, 2, "hello");
+  m_engine.seed(seed);
+  auto dist = std::normal_distribution<double>(0, 1);
   
   int nCreatedHits=0;
   int nDismissedHits=0;
@@ -186,7 +184,7 @@ DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits)
       // For this and other rng calls, use the fact that drawing from a
       // Gaussing with mean mu and variance sigma^2 is the same as drawing from
       // a normal distribution and then multiplying by sigma and adding mu
-      double tSmear  = resT > 0 ? rng() * resT : 0;
+      double tSmear  = resT > 0 ? dist(m_engine) * resT : 0;
       m_histograms[hT]->Fill(resT > 0 ? tSmear / resT : 0);
       m_histograms[diffT]->Fill(tSmear);
 
@@ -233,8 +231,8 @@ DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits)
       
       // if( tries > 0 ) debug() << "retry smearing for " <<  cellid_decoder( hit ).valueString() << " : retries " << tries << endmsg;
       
-      double uSmear  = rng() * resU;
-      double vSmear  = rng() * resV;
+      double uSmear  = dist(m_engine) * resU;
+      double vSmear  = dist(m_engine) * resV;
       
       dd4hep::rec::Vector3D newPosTmp;
       if (m_isStrip){
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
index d907bb5..ae0ed29 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
@@ -33,6 +33,7 @@
 #include "DDRec/SurfaceManager.h"
 
 #include <string>
+#include <random>
 #include <TH1F.h>
 
 /** ======= DDPlanarDigiProcessor ========== <br>
@@ -101,7 +102,7 @@ struct DDPlanarDigiProcessor final
   Gaudi::Property<std::string> m_geoSvcName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"};
   Gaudi::Property<int> m_maxTries{this, "MaxTries", 10, "Maximum number of tries to find a valid surface for a hit"};
 
-  Rndm::Numbers rng;
+  inline static thread_local std::mt19937 m_engine;
   const dd4hep::rec::SurfaceMap* surfaceMap;
   std::vector<TH1F*> m_histograms;
   SmartIF<IGeoSvc> m_geoSvc;

From f2cfe23beb249c5d8107920d29c7b2f4e5696551 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Tue, 3 Oct 2023 10:10:46 +0200
Subject: [PATCH 05/24] Add a m_collName to save the name of the collection

---
 .../components/DDPlanarDigiProcessor.cpp      | 21 +++++++++----------
 .../components/DDPlanarDigiProcessor.h        |  4 +++-
 2 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
index c66b622..4b09927 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
@@ -97,6 +97,9 @@ DDPlanarDigiProcessor::initialize() {
                                 << det.name() << " in SurfaceManager ";
     throw std::runtime_error(err.str());
   }
+
+  (void)this->getProperty("SimTrackerHitCollectionName", m_collName);
+
   return StatusCode::SUCCESS;
 }
 
@@ -104,11 +107,12 @@ DDPlanarDigiProcessor::initialize() {
 std::tuple<TrackerHitPlaneColl, Association>
 DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits) const {
 
-  // TODO: Set seed
-  // gslrng_set( rng, Global::EVENTSEEDER->getSeed(this) );   
-  // debug() << "seed set to " << Global::EVENTSEEDER->getSeed(this) << endmsg;
-  auto seed = m_uidSvc->getUniqueID(1, 2, "hello");
+  // TODO: Set seed correctly
+  auto seed = m_uidSvc->getUniqueID(1, 2, m_collName);
   m_engine.seed(seed);
+  // For rng calls, use the fact that drawing from a
+  // Gaussing with mean mu and variance sigma^2 is the same as drawing from
+  // a normal distribution and then multiplying by sigma and adding mu
   auto dist = std::normal_distribution<double>(0, 1);
   
   int nCreatedHits=0;
@@ -121,9 +125,7 @@ DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits)
   dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(cellIDEncodingString); 
   
   int nSimHits = simTrackerHits.size();
-  std::string collName;
-  (void)this->getProperty("SimTrackerHitCollectionName", collName);
-  debug() << "Processing collection " << collName << " with " << simTrackerHits.size() << " hits ... " << endmsg;
+  debug() << "Processing collection " << m_collName << " with " << simTrackerHits.size() << " hits ... " << endmsg;
 
   for(const auto& hit : simTrackerHits) {
     m_histograms[hitE]->Fill( hit.getEDep() * (dd4hep::GeV / dd4hep::keV) );
@@ -181,9 +183,6 @@ DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits)
     if (m_resTLayer.size() and m_resTLayer[0] > 0) {
       float resT = m_resTLayer.size() > 1 ? m_resTLayer[layer] : m_resTLayer[0];
 
-      // For this and other rng calls, use the fact that drawing from a
-      // Gaussing with mean mu and variance sigma^2 is the same as drawing from
-      // a normal distribution and then multiplying by sigma and adding mu
       double tSmear  = resT > 0 ? dist(m_engine) * resT : 0;
       m_histograms[hT]->Fill(resT > 0 ? tSmear / resT : 0);
       m_histograms[diffT]->Fill(tSmear);
@@ -215,7 +214,7 @@ DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits)
     dd4hep::rec::Vector3D u = surf->u();
     dd4hep::rec::Vector3D v = surf->v();
 
-    // get local coordinates on surface
+    // Get local coordinates on surface
     dd4hep::rec::Vector2D lv = surf->globalToLocal(dd4hep::mm * oldPos);
     double uL = lv[0] / dd4hep::mm;
     double vL = lv[1] / dd4hep::mm;
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
index ae0ed29..c557f37 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
@@ -102,9 +102,11 @@ struct DDPlanarDigiProcessor final
   Gaudi::Property<std::string> m_geoSvcName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"};
   Gaudi::Property<int> m_maxTries{this, "MaxTries", 10, "Maximum number of tries to find a valid surface for a hit"};
 
-  inline static thread_local std::mt19937 m_engine;
   const dd4hep::rec::SurfaceMap* surfaceMap;
   std::vector<TH1F*> m_histograms;
+  std::string m_collName;
+
+  inline static thread_local std::mt19937 m_engine;
   SmartIF<IGeoSvc> m_geoSvc;
   SmartIF<IUniqueIDGenSvc> m_uidSvc;
   

From feb51c553f25dc27c5e2a5ccb0fbefef8de7f7c0 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Wed, 25 Oct 2023 14:56:52 +0200
Subject: [PATCH 06/24] Add thread-safe histogramming

---
 CMakeLists.txt                                |   1 +
 .../components/DDPlanarDigiProcessor.cpp      | 392 +++++++++---------
 .../components/DDPlanarDigiProcessor.h        |  71 ++--
 .../options/runDDPlanarDigiProcessor.py       |   2 +-
 4 files changed, 242 insertions(+), 224 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index dcca7d4..795d5b0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,6 +33,7 @@ find_package(k4FWCore REQUIRED)
 find_package(Gaudi REQUIRED)
 find_package(DD4hep REQUIRED)
 find_package(k4SimGeant4 REQUIRED)
+find_package(k4geo REQUIRED)
 
 include(GNUInstallDirs)
 
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
index 4b09927..d265d6a 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
@@ -18,83 +18,76 @@
  */
 #include "DDPlanarDigiProcessor.h"
 
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
 #include "edm4hep/SimTrackerHit.h"
 #include "edm4hep/TrackerHitPlaneCollection.h"
-#include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
 
-#include "DD4hep/Detector.h"
-#include "DD4hep/DD4hepUnits.h"
+#include "Gaudi/Accumulators/Histogram.h"
+#include "Gaudi/Histograming/Sink/Utils.h"
+
 #include "DD4hep/BitFieldCoder.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/Detector.h"
 
-#include <TMath.h>
 #include <TFile.h>
+#include <TMath.h>
+#include <TProfile.h>
 
 #include <cmath>
+#include <iostream>
 #include <random>
 #include <sstream>
-#include <iostream>
 
-enum {
-  hu = 0,
-  hv,
-  hT,
-  hitE,
-  hitsAccepted,
-  diffu,
-  diffv,
-  diffT,
-  hSize 
-};
+enum { hu = 0, hv, hT, hitE, hitsAccepted, diffu, diffv, diffT, hSize };
 
 DDPlanarDigiProcessor::DDPlanarDigiProcessor(const std::string& name, ISvcLocator* svcLoc)
-    : MultiTransformer(
-          name, svcLoc,
-          {
-           KeyValue("SimTrackerHitCollectionName", "SimTrackerHitss"),
-           },
-          {KeyValue("TrackerHitCollectionName", "VTXTrackerHits"),
-           KeyValue("SimTrkHitRelCollection", "VTXTrackerHitRelations")}) {
-
-
+    : MultiTransformer(name, svcLoc,
+                       {
+                           KeyValue("SimTrackerHitCollectionName", "SimTrackerHits"),
+                           KeyValue("HeaderName", "EventHeader"),
+                       },
+                       {KeyValue("TrackerHitCollectionName", "VTXTrackerHits"),
+                        KeyValue("SimTrkHitRelCollection", "VTXTrackerHitRelations")}) {
   m_uidSvc = service<IUniqueIDGenSvc>("UniqueIDGenSvc", true);
   if (!m_uidSvc) {
     error() << "Unable to get UniqueIDGenSvc" << endmsg;
   }
 
-  if(m_resULayer.size() != m_resVLayer.size()) {
-    error() << "DDPlanarDigiProcessor - Inconsistent number of resolutions given for U and V coordinate: " 
-       << "ResolutionU  :" <<   m_resULayer.size() << " != ResolutionV : " <<  m_resVLayer.size();
+  if (m_resULayer.size() != m_resVLayer.size()) {
+    error() << "DDPlanarDigiProcessor - Inconsistent number of resolutions given for U and V coordinate: "
+            << "ResolutionU  :" << m_resULayer.size() << " != ResolutionV : " << m_resVLayer.size();
 
     throw std::runtime_error("DDPlanarDigiProcessor: Inconsistent number of resolutions given for U and V coordinate");
   }
 
   m_histograms.resize(hSize);
 
-  m_histograms[hu] = new TH1F("hu", "smearing u", 50, -5., +5.);
-  m_histograms[hv] = new TH1F("hv", "smearing v", 50, -5., +5.);
-  m_histograms[hT] = new TH1F("hT", "smearing time", 50, -5., +5.);
+  m_histograms[hu].reset(new Gaudi::Accumulators::Histogram<1>{this, "hu", "smearing u", {50, -5., +5.}});
+  m_histograms[hv].reset(new Gaudi::Accumulators::Histogram<1>{this, "hv", "smearing v", {50, -5., +5.}});
+  m_histograms[hT].reset(new Gaudi::Accumulators::Histogram<1>{this, "hT", "smearing time", {50, -5., +5.}});
 
-  m_histograms[diffu] = new TH1F("diffu", "diff u", 1000, -5., +5.);
-  m_histograms[diffv] = new TH1F("diffv", "diff v", 1000, -5., +5.);
-  m_histograms[diffT] = new TH1F("diffT", "diff time", 1000, -5., +5.);
+  m_histograms[diffu].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffu", "diff u", {1000, -.1, +.1}});
+  m_histograms[diffv].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffv", "diff v", {1000, -5., +5.}});
+  m_histograms[diffT].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffT", "diff time", {1000, -5., +5.}});
 
-  m_histograms[hitE] = new TH1F("hitE", "hitEnergy in keV", 1000, 0, 200);
-  m_histograms[hitsAccepted] = new TH1F("hitsAccepted", "Fraction of accepted hits [%]", 201, 0, 100.5);
+  m_histograms[hitE].reset(new Gaudi::Accumulators::Histogram<1>{this, "hitE", "hitEnergy in keV", {1000, 0, 200}});
+  m_histograms[hitsAccepted].reset(
+      new Gaudi::Accumulators::Histogram<1>{this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
 
   m_geoSvc = serviceLocator()->service(m_geoSvcName);
 }
 
-StatusCode
-DDPlanarDigiProcessor::initialize() {
+StatusCode DDPlanarDigiProcessor::initialize() {
   const auto detector = m_geoSvc->getDetector();
 
-  const auto surfMan = detector->extension<dd4hep::rec::SurfaceManager>();
-  dd4hep::DetElement det = detector->detector(m_subDetName.value());
-  surfaceMap = surfMan->map(m_subDetName.value());
+  const auto         surfMan = detector->extension<dd4hep::rec::SurfaceManager>();
+  dd4hep::DetElement det     = detector->detector(m_subDetName.value());
+  surfaceMap                 = surfMan->map(m_subDetName.value());
 
-  if(!surfaceMap) {   
-    std::stringstream err ; err << " Could not find surface map for detector: " 
-                                << det.name() << " in SurfaceManager ";
+  if (!surfaceMap) {
+    std::stringstream err;
+    err << " Could not find surface map for detector: " << det.name() << " in SurfaceManager ";
     throw std::runtime_error(err.str());
   }
 
@@ -103,34 +96,33 @@ DDPlanarDigiProcessor::initialize() {
   return StatusCode::SUCCESS;
 }
 
-
-std::tuple<TrackerHitPlaneColl, Association>
-DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits) const {
-
-  // TODO: Set seed correctly
-  auto seed = m_uidSvc->getUniqueID(1, 2, m_collName);
+std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigiProcessor::operator()(
+    const SimTrackerHitCollection& simTrackerHits, const Header& headers) const {
+  auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), m_collName);
+  info() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
+         << headers[0].getRunNumber() << endmsg;
   m_engine.seed(seed);
   // For rng calls, use the fact that drawing from a
-  // Gaussing with mean mu and variance sigma^2 is the same as drawing from
+  // Gaussian with mean mu and variance sigma^2 is the same as drawing from
   // a normal distribution and then multiplying by sigma and adding mu
   auto dist = std::normal_distribution<double>(0, 1);
-  
-  int nCreatedHits=0;
-  int nDismissedHits=0;
+
+  int nCreatedHits   = 0;
+  int nDismissedHits = 0;
 
   auto trkhitVec = edm4hep::TrackerHitPlaneCollection();
-  auto thsthcol = edm4hep::MCRecoTrackerHitPlaneAssociationCollection();
+  auto thsthcol  = edm4hep::MCRecoTrackerHitPlaneAssociationCollection();
 
   std::string cellIDEncodingString = m_geoSvc->constantAsString(m_encodingStringVariable.value());
-  dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(cellIDEncodingString); 
-  
+  dd4hep::DDSegmentation::BitFieldCoder bitFieldCoder(cellIDEncodingString);
+
   int nSimHits = simTrackerHits.size();
   debug() << "Processing collection " << m_collName << " with " << simTrackerHits.size() << " hits ... " << endmsg;
 
-  for(const auto& hit : simTrackerHits) {
-    m_histograms[hitE]->Fill( hit.getEDep() * (dd4hep::GeV / dd4hep::keV) );
+  for (const auto& hit : simTrackerHits) {
+    ++(*m_histograms[hitE])[hit.getEDep() * (dd4hep::GeV / dd4hep::keV)];
 
-    if(hit.getEDep() < m_minEnergy) {
+    if (hit.getEDep() < m_minEnergy) {
       debug() << "Hit with insufficient energy " << hit.getEDep() * (dd4hep::GeV / dd4hep::keV) << " keV" << endmsg;
       continue;
     }
@@ -140,215 +132,221 @@ DDPlanarDigiProcessor::operator()(const SimTrackerHitCollection& simTrackerHits)
     // get the measurement surface for this hit using the CellID
     dd4hep::rec::SurfaceMap::const_iterator sI = surfaceMap->find(cellID0);
 
-    if(sI == surfaceMap->end()){
-      std::stringstream err; err << " DDPlanarDigiProcessor::processEvent(): no surface found for cellID : " << std::endl;
-                                  // <<   cellid_decoder( hit ).valueString() ;
+    if (sI == surfaceMap->end()) {
+      std::stringstream err;
+      err << " DDPlanarDigiProcessor::processEvent(): no surface found for cellID : " << std::endl;
+      // <<   cellid_decoder( hit ).valueString() ;
       throw std::runtime_error(err.str());
     }
 
-    const dd4hep::rec::ISurface* surf = sI->second;
-    int layer = bitFieldCoder.get(cellID0, "layer");
+    const dd4hep::rec::ISurface* surf  = sI->second;
+    int                          layer = bitFieldCoder.get(cellID0, "layer");
 
     dd4hep::rec::Vector3D oldPos(hit.getPosition()[0], hit.getPosition()[1], hit.getPosition()[2]);
     dd4hep::rec::Vector3D newPos;
 
-    //  Check if Hit is inside sensitive 
+    //  Check if Hit is inside sensitive
     if (!surf->insideBounds(dd4hep::mm * oldPos)) {
-      
-      // debug() << "  hit at " << oldPos 
-      //         << " " << cellid_decoder( hit).valueString() 
-      //         << " is not on surface " 
-      //         << *surf  
+      // debug() << "  hit at " << oldPos
+      //         << " " << cellid_decoder( hit).valueString()
+      //         << " is not on surface "
+      //         << *surf
       //         << " distance: " << surf->distance(  dd4hep::mm * oldPos )
-      //         << endmsg;        
-
-      if(m_forceHitsOntoSurface){
-        dd4hep::rec::Vector2D lv = surf->globalToLocal(dd4hep::mm * oldPos);
-        dd4hep::rec::Vector3D oldPosOnSurf = (1./dd4hep::mm) * surf->localToGlobal(lv); 
-        
-        debug() << " moved to " << oldPosOnSurf << " distance " << (oldPosOnSurf-oldPos).r()
-                << endmsg;        
-          
+      //         << endmsg;
+
+      if (m_forceHitsOntoSurface) {
+        dd4hep::rec::Vector2D lv           = surf->globalToLocal(dd4hep::mm * oldPos);
+        dd4hep::rec::Vector3D oldPosOnSurf = (1. / dd4hep::mm) * surf->localToGlobal(lv);
+
+        debug() << " moved to " << oldPosOnSurf << " distance " << (oldPosOnSurf - oldPos).r() << endmsg;
+
         oldPos = oldPosOnSurf;
 
       } else {
         ++nDismissedHits;
-        continue; 
+        continue;
       }
     }
 
     // Smear time of the hit and apply the time window cut if needed
     double hitT = hit.getTime();
-    
+
     if (m_resTLayer.size() and m_resTLayer[0] > 0) {
       float resT = m_resTLayer.size() > 1 ? m_resTLayer[layer] : m_resTLayer[0];
 
-      double tSmear  = resT > 0 ? dist(m_engine) * resT : 0;
-      m_histograms[hT]->Fill(resT > 0 ? tSmear / resT : 0);
-      m_histograms[diffT]->Fill(tSmear);
+      double tSmear = resT > 0 ? dist(m_engine) * resT : 0;
+      ++(*m_histograms[hT])[resT > 0 ? tSmear / resT : 0];
+      ++(*m_histograms[diffT])[tSmear];
 
       hitT += tSmear;
-      debug() << "smeared hit at T: " << hit.getTime() << " ns to T: " << hitT << " ns according to resolution: " << resT << " ns" << endmsg;
+      debug() << "smeared hit at T: " << hit.getTime() << " ns to T: " << hitT
+              << " ns according to resolution: " << resT << " ns" << endmsg;
     }
-    
+
     // Correcting for the propagation time
     if (m_correctTimesForPropagation) {
-      double dt = oldPos.r() / (TMath::C()/1e6);
+      double dt = oldPos.r() / (TMath::C() / 1e6);
       hitT -= dt;
-      debug() << "corrected hit at R: " << oldPos.r() << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << endmsg;
+      debug() << "corrected hit at R: " << oldPos.r() << " mm by propagation time: " << dt << " ns to T: " << hitT
+              << " ns" << endmsg;
     }
-    
+
     // Skip the hit if its time is outside the acceptance time window
     if (m_useTimeWindow) {
       // TODO: Check that the length of the time window is OK
       float timeWindow_min = m_timeWindowMin.size() > 1 ? m_timeWindowMin[layer] : m_timeWindowMin[0];
       float timeWindow_max = m_timeWindowMax.size() > 1 ? m_timeWindowMax[layer] : m_timeWindowMax[0];
-      if ( hitT < timeWindow_min || hitT > timeWindow_max ) {
-        debug() << "hit at T: " << hit.getTime() << " smeared to: " << hitT << " is outside the time window: hit dropped"  << endmsg;
+      if (hitT < timeWindow_min || hitT > timeWindow_max) {
+        debug() << "hit at T: " << hit.getTime() << " smeared to: " << hitT
+                << " is outside the time window: hit dropped" << endmsg;
         ++nDismissedHits;
         continue;
       }
     }
-    
+
     // Try to smear the hit position but ensure the hit is inside the sensitive region
     dd4hep::rec::Vector3D u = surf->u();
     dd4hep::rec::Vector3D v = surf->v();
 
     // Get local coordinates on surface
     dd4hep::rec::Vector2D lv = surf->globalToLocal(dd4hep::mm * oldPos);
-    double uL = lv[0] / dd4hep::mm;
-    double vL = lv[1] / dd4hep::mm;
+    double                uL = lv[0] / dd4hep::mm;
+    double                vL = lv[1] / dd4hep::mm;
 
     bool acceptHit = false;
-    int tries =  0;              
-    
+    int  tries     = 0;
+
     // TODO: check lengths
     float resU = m_resULayer.size() > 1 ? m_resULayer[layer] : m_resULayer[0];
-    float resV = m_resVLayer.size() > 1 ? m_resVLayer[layer] : m_resVLayer[0]; 
+    float resV = m_resVLayer.size() > 1 ? m_resVLayer[layer] : m_resVLayer[0];
 
-    while(tries < m_maxTries) {
-      
+
+    while (tries < m_maxTries) {
       // if( tries > 0 ) debug() << "retry smearing for " <<  cellid_decoder( hit ).valueString() << " : retries " << tries << endmsg;
-      
-      double uSmear  = dist(m_engine) * resU;
-      double vSmear  = dist(m_engine) * resV;
-      
+
+      double uSmear = dist(m_engine) * resU;
+      double vSmear = dist(m_engine) * resV;
+
       dd4hep::rec::Vector3D newPosTmp;
-      if (m_isStrip){
-          if (m_subDetName == "SET"){
-            double xStripPos, yStripPos, zStripPos;
-            //Find intersection of the strip with the z=centerOfSensor plane to set it as the center of the SET strip
-            dd4hep::rec::Vector3D simHitPosSmeared = (1./dd4hep::mm) * ( surf->localToGlobal( dd4hep::rec::Vector2D( (uL+uSmear)*dd4hep::mm, 0.) ) );
-            zStripPos = surf->origin()[2] / dd4hep::mm;
-            double lineParam = (zStripPos - simHitPosSmeared[2])/v[2];
-            xStripPos = simHitPosSmeared[0] + lineParam*v[0];
-            yStripPos = simHitPosSmeared[1] + lineParam*v[1];
-            newPosTmp = dd4hep::rec::Vector3D(xStripPos, yStripPos, zStripPos);    
-          }
-          else{
-            newPosTmp = (1./dd4hep::mm) * ( surf->localToGlobal( dd4hep::rec::Vector2D( (uL+uSmear)*dd4hep::mm, 0. ) ) );    
-          }
-      }
-      else{
-        newPosTmp = (1./dd4hep::mm) * ( surf->localToGlobal( dd4hep::rec::Vector2D( (uL+uSmear)*dd4hep::mm, (vL+vSmear)*dd4hep::mm ) ) );
+      if (m_isStrip) {
+        if (m_subDetName == "SET") {
+          double xStripPos, yStripPos, zStripPos;
+          //Find intersection of the strip with the z=centerOfSensor plane to set it as the center of the SET strip
+          dd4hep::rec::Vector3D simHitPosSmeared =
+              (1. / dd4hep::mm) * (surf->localToGlobal(dd4hep::rec::Vector2D((uL + uSmear) * dd4hep::mm, 0.)));
+          zStripPos        = surf->origin()[2] / dd4hep::mm;
+          double lineParam = (zStripPos - simHitPosSmeared[2]) / v[2];
+          xStripPos        = simHitPosSmeared[0] + lineParam * v[0];
+          yStripPos        = simHitPosSmeared[1] + lineParam * v[1];
+          newPosTmp        = dd4hep::rec::Vector3D(xStripPos, yStripPos, zStripPos);
+        } else {
+          newPosTmp = (1. / dd4hep::mm) * (surf->localToGlobal(dd4hep::rec::Vector2D((uL + uSmear) * dd4hep::mm, 0.)));
+        }
+      } else {
+        newPosTmp =
+            (1. / dd4hep::mm) *
+            (surf->localToGlobal(dd4hep::rec::Vector2D((uL + uSmear) * dd4hep::mm, (vL + vSmear) * dd4hep::mm)));
       }
 
-      debug() << " hit at    : " << oldPos 
-              << " smeared to: " << newPosTmp
-              << " uL: " << uL 
-              << " vL: " << vL 
-              << " uSmear: " << uSmear
-              << " vSmear: " << vSmear
-              << endmsg;
-
+      debug() << " hit at    : " << oldPos << " smeared to: " << newPosTmp << " uL: " << uL << " vL: " << vL
+              << " uSmear: " << uSmear << " vSmear: " << vSmear << endmsg;
 
-      if ( surf->insideBounds( dd4hep::mm * newPosTmp ) ) {    
+      if (surf->insideBounds(dd4hep::mm * newPosTmp)) {
         acceptHit = true;
-        newPos     = newPosTmp;
+        newPos    = newPosTmp;
 
-        m_histograms[hu]->Fill(  uSmear / resU ); 
-        m_histograms[hv]->Fill(  vSmear / resV ); 
+        ++(*m_histograms[hu])[uSmear / resU];
+        ++(*m_histograms[hv])[vSmear / resV];
 
-        m_histograms[diffu]->Fill( uSmear );
-        m_histograms[diffv]->Fill( vSmear );
+        ++(*m_histograms[diffu])[uSmear];
+        ++(*m_histograms[diffv])[vSmear];
 
-        break;  
+        break;
       }
 
-      // debug() << "  hit at " << newPosTmp 
-      //         << " " << cellid_decoder( hit).valueString() 
-      //         << " is not on surface " 
-      //         << " distance: " << surf->distance( dd4hep::mm * newPosTmp ) 
-      //         << endmsg;        
-      
+      // debug() << "  hit at " << newPosTmp
+      //         << " " << cellid_decoder( hit).valueString()
+      //         << " is not on surface "
+      //         << " distance: " << surf->distance( dd4hep::mm * newPosTmp )
+      //         << endmsg;
+
       ++tries;
     }
 
-   if(!acceptHit) {
-     debug() << "hit could not be smeared within ladder after " << m_maxTries << "  tries: hit dropped"  << endmsg;
-     ++nDismissedHits;
-     continue; 
-   } 
-   
-   auto trkHit = trkhitVec.create();
-               
-   trkHit.setCellID( cellID0 );
-   
-   trkHit.setPosition( newPos.const_array()  );
-   trkHit.setTime( hitT );
-   trkHit.setEDep( hit.getEDep() );
-
-   float u_direction[2];
-   u_direction[0] = u.theta();
-   u_direction[1] = u.phi();
-   
-   float v_direction[2];
-   v_direction[0] = v.theta();
-   v_direction[1] = v.phi();
-
-   debug()  << " U[0] = "<< u_direction[0] << " U[1] = "<< u_direction[1] 
-                          << " V[0] = "<< v_direction[0] << " V[1] = "<< v_direction[1]
-                          << endmsg;
-
-   trkHit.setU(u_direction);
-   trkHit.setV(v_direction);
-   trkHit.setDu(resU);
-
-   if(m_isStrip) {
-     // store the resolution from the length of the wafer - in case a fitter might want to treat this as 2d hit ....
-     double stripRes = surf->length_along_v() / dd4hep::mm / std::sqrt(12);
-     trkHit.setDv(stripRes); 
-     // TODO: Set type?
-     // trkHit.setType( UTIL::set_bit( trkHit.getType() ,  UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ) );
-
-   } else {
-     trkHit.setDv(resV);
-   }
-
-   auto association = thsthcol.create();
-   association.setSim(hit);
-   association.setRec(trkHit);
-   
-   ++nCreatedHits;
-   
+    if (!acceptHit) {
+      debug() << "hit could not be smeared within ladder after " << m_maxTries << "  tries: hit dropped" << endmsg;
+      ++nDismissedHits;
+      continue;
+    }
+
+    auto trkHit = trkhitVec.create();
+
+    trkHit.setCellID(cellID0);
+
+    trkHit.setPosition(newPos.const_array());
+    trkHit.setTime(hitT);
+    trkHit.setEDep(hit.getEDep());
+
+    float u_direction[2];
+    u_direction[0] = u.theta();
+    u_direction[1] = u.phi();
+
+    float v_direction[2];
+    v_direction[0] = v.theta();
+    v_direction[1] = v.phi();
+
+    debug() << " U[0] = " << u_direction[0] << " U[1] = " << u_direction[1] << " V[0] = " << v_direction[0]
+            << " V[1] = " << v_direction[1] << endmsg;
+
+    trkHit.setU(u_direction);
+    trkHit.setV(v_direction);
+    trkHit.setDu(resU);
+
+    if (m_isStrip) {
+      // store the resolution from the length of the wafer - in case a fitter might want to treat this as 2d hit ....
+      double stripRes = surf->length_along_v() / dd4hep::mm / std::sqrt(12);
+      trkHit.setDv(stripRes);
+      // TODO: Set type?
+      // trkHit.setType( UTIL::set_bit( trkHit.getType() ,  UTIL::ILDTrkHitTypeBit::ONE_DIMENSIONAL ) );
+
+    } else {
+      trkHit.setDv(resV);
+    }
+
+    auto association = thsthcol.create();
+    association.setSim(hit);
+    association.setRec(trkHit);
+
+    ++nCreatedHits;
   }
 
   // Filling the fraction of accepted hits in the event
   float accFraction = nSimHits > 0 ? float(nCreatedHits) / float(nSimHits) * 100.0 : 0.0;
-  m_histograms[hitsAccepted]->Fill(accFraction);
-  
+  ++(*m_histograms[hitsAccepted])[accFraction];
+
   debug() << "Created " << nCreatedHits << " hits, " << nDismissedHits << " hits  dismissed" << endmsg;
-  
+
   return std::make_tuple(std::move(trkhitVec), std::move(thsthcol));
 }
 
-
-StatusCode
-DDPlanarDigiProcessor::finalize() {
+StatusCode DDPlanarDigiProcessor::finalize() {
   auto file = TFile::Open("DDPlanarDigiProcessor.root", "RECREATE");
+  auto names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
+  auto it = names.begin();
   for (auto& h : m_histograms) {
-    h->Write();
-    delete h;
+    std::string name = "";
+    // Name that will appear in the stats table
+    std::string    histName = *it;
+    nlohmann::json json = *h;
+    // Name of the histogram in the ROOT file
+    if (strcmp(*it, "diffu") == 0) {
+      info() << json << endmsg;
+    }
+    auto [histo, dir] =
+        Gaudi::Histograming::Sink::jsonToRootHistogram<Gaudi::Histograming::Sink::Traits<false, TH1D, 1>>(
+            name, histName, json);
+    histo.Write(*it);
+    ++it;
   }
   file->Close();
   return StatusCode::SUCCESS;
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
index c557f37..f4f3cfa 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
@@ -19,22 +19,24 @@
 #ifndef DDPLANARDIGIPROCESSOR_H
 #define DDPLANARDIGIPROCESSOR_H
 
+#include "Gaudi/Accumulators/Histogram.h"
 #include "Gaudi/Property.h"
 #include "GaudiAlg/Transformer.h"
 #include "GaudiKernel/RndmGenerators.h"
-#include "k4Interface/IGeoSvc.h"
 #include "k4FWCore/BaseClass.h"
+#include "k4Interface/IGeoSvc.h"
 #include "k4Interface/IUniqueIDGenSvc.h"
 
+#include "edm4hep/EventHeaderCollection.h"
+#include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/TrackerHitPlaneCollection.h"
-#include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
 
 #include "DDRec/SurfaceManager.h"
 
-#include <string>
-#include <random>
 #include <TH1F.h>
+#include <random>
+#include <string>
 
 /** ======= DDPlanarDigiProcessor ========== <br>
  * Creates TrackerHits from SimTrackerHits, smearing them according to the input parameters. 
@@ -69,33 +71,49 @@
  */
 
 using SimTrackerHitCollection = edm4hep::SimTrackerHitCollection;
+using Header                  = edm4hep::EventHeaderCollection;
 
 using TrackerHitPlaneColl = edm4hep::TrackerHitPlaneCollection;
-using Association = edm4hep::MCRecoTrackerHitPlaneAssociationCollection;
+using Association         = edm4hep::MCRecoTrackerHitPlaneAssociationCollection;
 
 struct DDPlanarDigiProcessor final
-    : Gaudi::Functional::MultiTransformer<std::tuple<TrackerHitPlaneColl, Association>(const SimTrackerHitCollection&), BaseClass_t> {
-
+    : Gaudi::Functional::MultiTransformer<
+          std::tuple<TrackerHitPlaneColl, Association>(const SimTrackerHitCollection&, const Header&), BaseClass_t> {
   DDPlanarDigiProcessor(const std::string& name, ISvcLocator* svcLoc);
-  
+
   StatusCode initialize() override;
   StatusCode finalize() override;
 
-  std::tuple<TrackerHitPlaneColl, Association> operator()(const SimTrackerHitCollection& simTrackerHits) const override;
-
- private:
-  Gaudi::Property<std::string> m_subDetName{this, "SubDetectorName", "VXD", "Name of the subdetector"};
-  Gaudi::Property<bool> m_isStrip{this, "IsStrip", false, "Whether the hits are 1D strip hits"};
-  Gaudi::Property<std::vector<float>> m_resULayer{this, "ResolutionU", {0.004}, "Resolution in the direction of u; either one per layer or one for all layers"};
-  Gaudi::Property<std::vector<float>> m_resVLayer{this, "ResolutionV", {0.004}, "Resolution in the direction of v; either one per layer or one for all layers"};
-  Gaudi::Property<std::vector<float>> m_resTLayer{this, "ResolutionT", {0.004}, "Resolution in the direction of t; either one per layer or one for all layers. If the single entry is negative, disable time smearing. "};
-  Gaudi::Property<bool> m_forceHitsOntoSurface{this, "ForceHitsOntoSurface", false, "Project hits onto the surfoce in case they are not yet on the surface"};
+  std::tuple<TrackerHitPlaneColl, Association> operator()(const SimTrackerHitCollection& simTrackerHits,
+                                                          const Header&                  headers) const override;
+
+private:
+  Gaudi::Property<std::string>        m_subDetName{this, "SubDetectorName", "VXD", "Name of the subdetector"};
+  Gaudi::Property<bool>               m_isStrip{this, "IsStrip", false, "Whether the hits are 1D strip hits"};
+  Gaudi::Property<std::vector<float>> m_resULayer{
+      this, "ResolutionU", {0.004}, "Resolution in the direction of u; either one per layer or one for all layers"};
+  Gaudi::Property<std::vector<float>> m_resVLayer{
+      this, "ResolutionV", {0.004}, "Resolution in the direction of v; either one per layer or one for all layers"};
+  Gaudi::Property<std::vector<float>> m_resTLayer{
+      this,
+      "ResolutionT",
+      {0.004},
+      "Resolution in the direction of t; either one per layer or one for all layers. If the single entry is negative, "
+      "disable time smearing. "};
+  Gaudi::Property<bool>   m_forceHitsOntoSurface{this, "ForceHitsOntoSurface", false,
+                                               "Project hits onto the surfoce in case they are not yet on the surface"};
   Gaudi::Property<double> m_minEnergy{this, "MinEnergy", 0.0, "Minimum energy (GeV) of SimTrackerHit to be digitized"};
 
-  Gaudi::Property<bool> m_useTimeWindow{this, "UseTimeWindow", false, "Only accept hits with time (after smearing) within the specified time window (default: false)"};
-  Gaudi::Property<bool> m_correctTimesForPropagation{this, "CorrectTimesForPropagation", false, "Correct hit time for the propagation: radial distance/c (default: false)"};
-  Gaudi::Property<std::vector<float>> m_timeWindowMin{this, "TimeWindowMin", {-1e9}, "Minimum time (ns) of SimTrackerHit to be digitized"};
-  Gaudi::Property<std::vector<float>> m_timeWindowMax{this, "TimeWindowMax", {1e9}, "Maximum time (ns) of SimTrackerHit to be digitized"};
+  Gaudi::Property<bool> m_useTimeWindow{
+      this, "UseTimeWindow", false,
+      "Only accept hits with time (after smearing) within the specified time window (default: false)"};
+  Gaudi::Property<bool> m_correctTimesForPropagation{
+      this, "CorrectTimesForPropagation", false,
+      "Correct hit time for the propagation: radial distance/c (default: false)"};
+  Gaudi::Property<std::vector<float>> m_timeWindowMin{
+      this, "TimeWindowMin", {-1e9}, "Minimum time (ns) of SimTrackerHit to be digitized"};
+  Gaudi::Property<std::vector<float>> m_timeWindowMax{
+      this, "TimeWindowMax", {1e9}, "Maximum time (ns) of SimTrackerHit to be digitized"};
   Gaudi::Property<std::string> m_encodingStringVariable{
       this, "EncodingStringParameterName", "GlobalTrackerReadoutID",
       "The name of the DD4hep constant that contains the Encoding string for tracking detectors"};
@@ -103,13 +121,14 @@ struct DDPlanarDigiProcessor final
   Gaudi::Property<int> m_maxTries{this, "MaxTries", 10, "Maximum number of tries to find a valid surface for a hit"};
 
   const dd4hep::rec::SurfaceMap* surfaceMap;
-  std::vector<TH1F*> m_histograms;
-  std::string m_collName;
+  // std::vector<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>> m_histograms;
+  std::vector<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>> m_histograms;
+  std::string                                                     m_collName;
 
+  // inline static thread_local std::mt19937 m_engine;
   inline static thread_local std::mt19937 m_engine;
-  SmartIF<IGeoSvc> m_geoSvc;
-  SmartIF<IUniqueIDGenSvc> m_uidSvc;
-  
+  SmartIF<IGeoSvc>                        m_geoSvc;
+  SmartIF<IUniqueIDGenSvc>                m_uidSvc;
 };
 
 DECLARE_COMPONENT(DDPlanarDigiProcessor)
diff --git a/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py b/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
index 8ac0250..07a4988 100644
--- a/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
+++ b/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
@@ -45,6 +45,7 @@
 inp = PodioInput()
 inp.collections = [
     "VertexBarrelCollection",
+    "EventHeader",
 ]
 
 out = PodioOutput("out")
@@ -56,4 +57,3 @@
                ExtSvc=[data_svc],
                OutputLevel=INFO,
                )
-

From 08aae5060febb62033f452bba1ec07ef0c5ce1ed Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Wed, 25 Oct 2023 15:01:45 +0200
Subject: [PATCH 07/24] Add some default options and remove some includes

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp | 4 +---
 k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h   | 4 +++-
 k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py  | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
index d265d6a..73e1e13 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
@@ -31,8 +31,6 @@
 #include "DD4hep/Detector.h"
 
 #include <TFile.h>
-#include <TMath.h>
-#include <TProfile.h>
 
 #include <cmath>
 #include <iostream>
@@ -330,7 +328,7 @@ std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigiProcessor::operator()(
 }
 
 StatusCode DDPlanarDigiProcessor::finalize() {
-  auto file = TFile::Open("DDPlanarDigiProcessor.root", "RECREATE");
+  auto file = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
   auto names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
   auto it = names.begin();
   for (auto& h : m_histograms) {
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
index f4f3cfa..8087fdb 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
@@ -34,9 +34,9 @@
 
 #include "DDRec/SurfaceManager.h"
 
-#include <TH1F.h>
 #include <random>
 #include <string>
+#include <vector>
 
 /** ======= DDPlanarDigiProcessor ========== <br>
  * Creates TrackerHits from SimTrackerHits, smearing them according to the input parameters. 
@@ -120,6 +120,8 @@ struct DDPlanarDigiProcessor final
   Gaudi::Property<std::string> m_geoSvcName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"};
   Gaudi::Property<int> m_maxTries{this, "MaxTries", 10, "Maximum number of tries to find a valid surface for a hit"};
 
+  Gaudi::Property<std::string> m_outputFileName{this, "OutputFileName", "planar_digi_histograms.root", "Output file name for the histograms"};
+
   const dd4hep::rec::SurfaceMap* surfaceMap;
   // std::vector<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>> m_histograms;
   std::vector<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>> m_histograms;
diff --git a/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py b/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
index 07a4988..a9dc1b1 100644
--- a/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
+++ b/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
@@ -49,7 +49,7 @@
 ]
 
 out = PodioOutput("out")
-out.filename = "ttbar_edm4hep_digi.root"
+out.filename = "planar_digi_histograms.root"
 
 ApplicationMgr(TopAlg=[inp, processor, out],
                EvtSel="NONE",

From 6a221caa93d4b379daaea5eebcbc01901cc1ae3d Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Wed, 25 Oct 2023 15:13:44 +0200
Subject: [PATCH 08/24] Remove processor from the name

---
 k4Reco/CMakeLists.txt                            |  2 +-
 ...DPlanarDigiProcessor.cpp => DDPlanarDigi.cpp} | 16 ++++++++--------
 .../{DDPlanarDigiProcessor.h => DDPlanarDigi.h}  | 14 +++++++-------
 ...PlanarDigiProcessor.py => runDDPlanarDigi.py} |  4 ++--
 4 files changed, 18 insertions(+), 18 deletions(-)
 rename k4Reco/DDPlanarDigi/components/{DDPlanarDigiProcessor.cpp => DDPlanarDigi.cpp} (95%)
 rename k4Reco/DDPlanarDigi/components/{DDPlanarDigiProcessor.h => DDPlanarDigi.h} (95%)
 rename k4Reco/DDPlanarDigi/options/{runDDPlanarDigiProcessor.py => runDDPlanarDigi.py} (95%)

diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt
index 35fd9c9..3c35720 100644
--- a/k4Reco/CMakeLists.txt
+++ b/k4Reco/CMakeLists.txt
@@ -17,7 +17,7 @@ See the License for the specific language governing permissions and
 limitations under the License.
 ]]
 
-set(sources DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+set(sources DDPlanarDigi/components/DDPlanarDigi.cpp
 )
 gaudi_add_module(k4RecoPlugins
                  SOURCES ${sources}
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
similarity index 95%
rename from k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
rename to k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 73e1e13..cc2fed0 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -16,7 +16,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-#include "DDPlanarDigiProcessor.h"
+#include "DDPlanarDigi.h"
 
 #include "edm4hep/EventHeaderCollection.h"
 #include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
@@ -39,7 +39,7 @@
 
 enum { hu = 0, hv, hT, hitE, hitsAccepted, diffu, diffv, diffT, hSize };
 
-DDPlanarDigiProcessor::DDPlanarDigiProcessor(const std::string& name, ISvcLocator* svcLoc)
+DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     : MultiTransformer(name, svcLoc,
                        {
                            KeyValue("SimTrackerHitCollectionName", "SimTrackerHits"),
@@ -53,10 +53,10 @@ DDPlanarDigiProcessor::DDPlanarDigiProcessor(const std::string& name, ISvcLocato
   }
 
   if (m_resULayer.size() != m_resVLayer.size()) {
-    error() << "DDPlanarDigiProcessor - Inconsistent number of resolutions given for U and V coordinate: "
+    error() << "DDPlanarDigi - Inconsistent number of resolutions given for U and V coordinate: "
             << "ResolutionU  :" << m_resULayer.size() << " != ResolutionV : " << m_resVLayer.size();
 
-    throw std::runtime_error("DDPlanarDigiProcessor: Inconsistent number of resolutions given for U and V coordinate");
+    throw std::runtime_error("DDPlanarDigi: Inconsistent number of resolutions given for U and V coordinate");
   }
 
   m_histograms.resize(hSize);
@@ -76,7 +76,7 @@ DDPlanarDigiProcessor::DDPlanarDigiProcessor(const std::string& name, ISvcLocato
   m_geoSvc = serviceLocator()->service(m_geoSvcName);
 }
 
-StatusCode DDPlanarDigiProcessor::initialize() {
+StatusCode DDPlanarDigi::initialize() {
   const auto detector = m_geoSvc->getDetector();
 
   const auto         surfMan = detector->extension<dd4hep::rec::SurfaceManager>();
@@ -94,7 +94,7 @@ StatusCode DDPlanarDigiProcessor::initialize() {
   return StatusCode::SUCCESS;
 }
 
-std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigiProcessor::operator()(
+std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigi::operator()(
     const SimTrackerHitCollection& simTrackerHits, const Header& headers) const {
   auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), m_collName);
   info() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
@@ -132,7 +132,7 @@ std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigiProcessor::operator()(
 
     if (sI == surfaceMap->end()) {
       std::stringstream err;
-      err << " DDPlanarDigiProcessor::processEvent(): no surface found for cellID : " << std::endl;
+      err << " DDPlanarDigi::processEvent(): no surface found for cellID : " << std::endl;
       // <<   cellid_decoder( hit ).valueString() ;
       throw std::runtime_error(err.str());
     }
@@ -327,7 +327,7 @@ std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigiProcessor::operator()(
   return std::make_tuple(std::move(trkhitVec), std::move(thsthcol));
 }
 
-StatusCode DDPlanarDigiProcessor::finalize() {
+StatusCode DDPlanarDigi::finalize() {
   auto file = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
   auto names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
   auto it = names.begin();
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
similarity index 95%
rename from k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
rename to k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 8087fdb..3e8b808 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigiProcessor.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -16,8 +16,8 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-#ifndef DDPLANARDIGIPROCESSOR_H
-#define DDPLANARDIGIPROCESSOR_H
+#ifndef DDPLANARDIGI_H
+#define DDPLANARDIGI_H
 
 #include "Gaudi/Accumulators/Histogram.h"
 #include "Gaudi/Property.h"
@@ -38,7 +38,7 @@
 #include <string>
 #include <vector>
 
-/** ======= DDPlanarDigiProcessor ========== <br>
+/** ======= DDPlanarDigi ========== <br>
  * Creates TrackerHits from SimTrackerHits, smearing them according to the input parameters. 
  * The positions of "digitized" TrackerHits are obtained by gaussian smearing positions
  * of SimTrackerHits perpendicular and along the ladder according to the specified point resolutions. 
@@ -67,7 +67,7 @@
  * 
  * @author F.Gaede CERN/DESY, S. Aplin DESY
  * @date Dec 2014
- * Originally in https://github.com/iLCSoft/MarlinTrkProcessors/blob/master/source/Digitisers/include/DDPlanarDigiProcessor.h
+ * Originally in https://github.com/iLCSoft/MarlinTrkProcessors/blob/master/source/Digitisers/include/DDPlanarDigi.h
  */
 
 using SimTrackerHitCollection = edm4hep::SimTrackerHitCollection;
@@ -76,10 +76,10 @@ using Header                  = edm4hep::EventHeaderCollection;
 using TrackerHitPlaneColl = edm4hep::TrackerHitPlaneCollection;
 using Association         = edm4hep::MCRecoTrackerHitPlaneAssociationCollection;
 
-struct DDPlanarDigiProcessor final
+struct DDPlanarDigi final
     : Gaudi::Functional::MultiTransformer<
           std::tuple<TrackerHitPlaneColl, Association>(const SimTrackerHitCollection&, const Header&), BaseClass_t> {
-  DDPlanarDigiProcessor(const std::string& name, ISvcLocator* svcLoc);
+  DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc);
 
   StatusCode initialize() override;
   StatusCode finalize() override;
@@ -133,6 +133,6 @@ struct DDPlanarDigiProcessor final
   SmartIF<IUniqueIDGenSvc>                m_uidSvc;
 };
 
-DECLARE_COMPONENT(DDPlanarDigiProcessor)
+DECLARE_COMPONENT(DDPlanarDigi)
 
 #endif
diff --git a/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py b/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
similarity index 95%
rename from k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
rename to k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
index a9dc1b1..16358eb 100644
--- a/k4Reco/DDPlanarDigi/options/runDDPlanarDigiProcessor.py
+++ b/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
@@ -18,7 +18,7 @@
 #
 from Gaudi.Configuration import INFO
 from Configurables import ApplicationMgr, k4DataSvc, PodioInput, PodioOutput
-from Configurables import DDPlanarDigiProcessor
+from Configurables import DDPlanarDigi
 from Configurables import GeoSvc
 from Configurables import UniqueIDGenSvc
 import os
@@ -30,7 +30,7 @@
 geoservice.OutputLevel = INFO
 geoservice.EnableGeant4Geo = False
 
-processor = DDPlanarDigiProcessor()
+processor = DDPlanarDigi()
 processor.SubDetectorName = "Vertex"
 processor.IsStrip = False
 processor.ResolutionU = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]

From d985a11ee1e57c0710c3a7066b28fcfbdd7b9446 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Wed, 25 Oct 2023 21:25:44 +0200
Subject: [PATCH 09/24] Change the bins for diffv

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index cc2fed0..d590f7f 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -66,7 +66,7 @@ DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
   m_histograms[hT].reset(new Gaudi::Accumulators::Histogram<1>{this, "hT", "smearing time", {50, -5., +5.}});
 
   m_histograms[diffu].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffu", "diff u", {1000, -.1, +.1}});
-  m_histograms[diffv].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffv", "diff v", {1000, -5., +5.}});
+  m_histograms[diffv].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffv", "diff v", {1000, -.1, +.1}});
   m_histograms[diffT].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffT", "diff time", {1000, -5., +5.}});
 
   m_histograms[hitE].reset(new Gaudi::Accumulators::Histogram<1>{this, "hitE", "hitEnergy in keV", {1000, 0, 200}});

From ee0469434b0e57e5e07da0abbe0e944244c7f486 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Thu, 26 Oct 2023 09:54:00 +0200
Subject: [PATCH 10/24] Remove info message

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index d590f7f..258e887 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -336,10 +336,6 @@ StatusCode DDPlanarDigi::finalize() {
     // Name that will appear in the stats table
     std::string    histName = *it;
     nlohmann::json json = *h;
-    // Name of the histogram in the ROOT file
-    if (strcmp(*it, "diffu") == 0) {
-      info() << json << endmsg;
-    }
     auto [histo, dir] =
         Gaudi::Histograming::Sink::jsonToRootHistogram<Gaudi::Histograming::Sink::Traits<false, TH1D, 1>>(
             name, histName, json);

From feb1044607efa0c3277772f8099ec581026933c1 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Fri, 27 Oct 2023 10:55:32 +0200
Subject: [PATCH 11/24] Use the name of the algorithm

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 258e887..ea5fe5d 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -89,6 +89,7 @@ StatusCode DDPlanarDigi::initialize() {
     throw std::runtime_error(err.str());
   }
 
+  // Get and store the name for a debug message later
   (void)this->getProperty("SimTrackerHitCollectionName", m_collName);
 
   return StatusCode::SUCCESS;
@@ -96,7 +97,7 @@ StatusCode DDPlanarDigi::initialize() {
 
 std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigi::operator()(
     const SimTrackerHitCollection& simTrackerHits, const Header& headers) const {
-  auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), m_collName);
+  auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), this->name());
   info() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
          << headers[0].getRunNumber() << endmsg;
   m_engine.seed(seed);

From 4351e68bbc525918fe21119f7b9eb13d4df9a538 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Mon, 30 Oct 2023 18:38:55 +0100
Subject: [PATCH 12/24] Address PR comments

---
 CMakeLists.txt                                |  1 -
 .../DDPlanarDigi/components/DDPlanarDigi.cpp  |  3 ---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h | 26 +++++++++----------
 test/CMakeLists.txt                           |  2 ++
 4 files changed, 14 insertions(+), 18 deletions(-)

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 795d5b0..dcca7d4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,7 +33,6 @@ find_package(k4FWCore REQUIRED)
 find_package(Gaudi REQUIRED)
 find_package(DD4hep REQUIRED)
 find_package(k4SimGeant4 REQUIRED)
-find_package(k4geo REQUIRED)
 
 include(GNUInstallDirs)
 
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index ea5fe5d..9a444bf 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -37,7 +37,6 @@
 #include <random>
 #include <sstream>
 
-enum { hu = 0, hv, hT, hitE, hitsAccepted, diffu, diffv, diffT, hSize };
 
 DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     : MultiTransformer(name, svcLoc,
@@ -59,8 +58,6 @@ DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     throw std::runtime_error("DDPlanarDigi: Inconsistent number of resolutions given for U and V coordinate");
   }
 
-  m_histograms.resize(hSize);
-
   m_histograms[hu].reset(new Gaudi::Accumulators::Histogram<1>{this, "hu", "smearing u", {50, -5., +5.}});
   m_histograms[hv].reset(new Gaudi::Accumulators::Histogram<1>{this, "hv", "smearing v", {50, -5., +5.}});
   m_histograms[hT].reset(new Gaudi::Accumulators::Histogram<1>{this, "hT", "smearing time", {50, -5., +5.}});
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 3e8b808..763fc38 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -70,22 +70,21 @@
  * Originally in https://github.com/iLCSoft/MarlinTrkProcessors/blob/master/source/Digitisers/include/DDPlanarDigi.h
  */
 
-using SimTrackerHitCollection = edm4hep::SimTrackerHitCollection;
-using Header                  = edm4hep::EventHeaderCollection;
-
-using TrackerHitPlaneColl = edm4hep::TrackerHitPlaneCollection;
-using Association         = edm4hep::MCRecoTrackerHitPlaneAssociationCollection;
+enum { hu = 0, hv, hT, hitE, hitsAccepted, diffu, diffv, diffT, hSize };
 
 struct DDPlanarDigi final
     : Gaudi::Functional::MultiTransformer<
-          std::tuple<TrackerHitPlaneColl, Association>(const SimTrackerHitCollection&, const Header&), BaseClass_t> {
+          std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAssociationCollection>(
+              const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&),
+          BaseClass_t> {
   DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc);
 
   StatusCode initialize() override;
   StatusCode finalize() override;
 
-  std::tuple<TrackerHitPlaneColl, Association> operator()(const SimTrackerHitCollection& simTrackerHits,
-                                                          const Header&                  headers) const override;
+  std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAssociationCollection> operator()(
+      const edm4hep::SimTrackerHitCollection& simTrackerHits,
+      const edm4hep::EventHeaderCollection&   headers) const override;
 
 private:
   Gaudi::Property<std::string>        m_subDetName{this, "SubDetectorName", "VXD", "Name of the subdetector"};
@@ -120,14 +119,13 @@ struct DDPlanarDigi final
   Gaudi::Property<std::string> m_geoSvcName{this, "GeoSvcName", "GeoSvc", "The name of the GeoSvc instance"};
   Gaudi::Property<int> m_maxTries{this, "MaxTries", 10, "Maximum number of tries to find a valid surface for a hit"};
 
-  Gaudi::Property<std::string> m_outputFileName{this, "OutputFileName", "planar_digi_histograms.root", "Output file name for the histograms"};
+  Gaudi::Property<std::string> m_outputFileName{this, "OutputFileName", "planar_digi_histograms.root",
+                                                "Output file name for the histograms"};
 
-  const dd4hep::rec::SurfaceMap* surfaceMap;
-  // std::vector<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>> m_histograms;
-  std::vector<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>> m_histograms;
-  std::string                                                     m_collName;
+  const dd4hep::rec::SurfaceMap*                                        surfaceMap;
+  std::array<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>, hSize> m_histograms;
+  std::string                                                           m_collName;
 
-  // inline static thread_local std::mt19937 m_engine;
   inline static thread_local std::mt19937 m_engine;
   SmartIF<IGeoSvc>                        m_geoSvc;
   SmartIF<IUniqueIDGenSvc>                m_uidSvc;
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index c0b1b51..a2d02c8 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -17,3 +17,5 @@ See the License for the specific language governing permissions and
 limitations under the License.
 ]]
 
+
+find_package(k4geo REQUIRED)

From 82e657dfc5740e7bea30a0793aadbf8a42785419 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Tue, 31 Oct 2023 09:30:59 +0100
Subject: [PATCH 13/24] Fix more aliases

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 9a444bf..6df422c 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -92,8 +92,8 @@ StatusCode DDPlanarDigi::initialize() {
   return StatusCode::SUCCESS;
 }
 
-std::tuple<TrackerHitPlaneColl, Association> DDPlanarDigi::operator()(
-    const SimTrackerHitCollection& simTrackerHits, const Header& headers) const {
+std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAssociationCollection> DDPlanarDigi::operator()(
+    const edm4hep::SimTrackerHitCollection& simTrackerHits, const edm4hep::EventHeaderCollection& headers) const {
   auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), this->name());
   info() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
          << headers[0].getRunNumber() << endmsg;

From 895576a1e809fe7938cddd99fe0389aff4d89a51 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Tue, 31 Oct 2023 14:42:23 +0100
Subject: [PATCH 14/24] Add a sigma counter to have accurate values of the mean
 and StDev

---
 .../DDPlanarDigi/components/DDPlanarDigi.cpp  | 37 +++++++++++++++----
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h |  6 ++-
 2 files changed, 34 insertions(+), 9 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 6df422c..f37ad56 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -37,7 +37,6 @@
 #include <random>
 #include <sstream>
 
-
 DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     : MultiTransformer(name, svcLoc,
                        {
@@ -70,6 +69,16 @@ DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
   m_histograms[hitsAccepted].reset(
       new Gaudi::Accumulators::Histogram<1>{this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
 
+  m_sigma[hu].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+  m_sigma[hv].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+  m_sigma[hT].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+  m_sigma[diffu].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+  m_sigma[diffv].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+  m_sigma[diffT].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+  m_sigma[hitE].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+  m_sigma[hitsAccepted].reset(
+      new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+
   m_geoSvc = serviceLocator()->service(m_geoSvcName);
 }
 
@@ -92,11 +101,12 @@ StatusCode DDPlanarDigi::initialize() {
   return StatusCode::SUCCESS;
 }
 
-std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAssociationCollection> DDPlanarDigi::operator()(
-    const edm4hep::SimTrackerHitCollection& simTrackerHits, const edm4hep::EventHeaderCollection& headers) const {
+std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAssociationCollection>
+DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
+                         const edm4hep::EventHeaderCollection&   headers) const {
   auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), this->name());
-  info() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
-         << headers[0].getRunNumber() << endmsg;
+  debug() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
+          << headers[0].getRunNumber() << endmsg;
   m_engine.seed(seed);
   // For rng calls, use the fact that drawing from a
   // Gaussian with mean mu and variance sigma^2 is the same as drawing from
@@ -117,6 +127,7 @@ std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAss
 
   for (const auto& hit : simTrackerHits) {
     ++(*m_histograms[hitE])[hit.getEDep() * (dd4hep::GeV / dd4hep::keV)];
+    *m_sigma[hitE] += hit.getEDep() * (dd4hep::GeV / dd4hep::keV);
 
     if (hit.getEDep() < m_minEnergy) {
       debug() << "Hit with insufficient energy " << hit.getEDep() * (dd4hep::GeV / dd4hep::keV) << " keV" << endmsg;
@@ -172,7 +183,9 @@ std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAss
 
       double tSmear = resT > 0 ? dist(m_engine) * resT : 0;
       ++(*m_histograms[hT])[resT > 0 ? tSmear / resT : 0];
+      *m_sigma[hT] += resT > 0 ? tSmear / resT : 0;
       ++(*m_histograms[diffT])[tSmear];
+      *m_sigma[diffT] += tSmear;
 
       hitT += tSmear;
       debug() << "smeared hit at T: " << hit.getTime() << " ns to T: " << hitT
@@ -216,7 +229,6 @@ std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAss
     float resU = m_resULayer.size() > 1 ? m_resULayer[layer] : m_resULayer[0];
     float resV = m_resVLayer.size() > 1 ? m_resVLayer[layer] : m_resVLayer[0];
 
-
     while (tries < m_maxTries) {
       // if( tries > 0 ) debug() << "retry smearing for " <<  cellid_decoder( hit ).valueString() << " : retries " << tries << endmsg;
 
@@ -252,10 +264,14 @@ std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAss
         newPos    = newPosTmp;
 
         ++(*m_histograms[hu])[uSmear / resU];
+        *m_sigma[hu] += uSmear / resU;
         ++(*m_histograms[hv])[vSmear / resV];
+        *m_sigma[hv] += vSmear / resV;
 
         ++(*m_histograms[diffu])[uSmear];
+        *m_sigma[diffu] += uSmear;
         ++(*m_histograms[diffv])[vSmear];
+        *m_sigma[diffv] += vSmear;
 
         break;
       }
@@ -319,6 +335,7 @@ std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAss
   // Filling the fraction of accepted hits in the event
   float accFraction = nSimHits > 0 ? float(nCreatedHits) / float(nSimHits) * 100.0 : 0.0;
   ++(*m_histograms[hitsAccepted])[accFraction];
+  *m_sigma[hitsAccepted] += accFraction;
 
   debug() << "Created " << nCreatedHits << " hits, " << nDismissedHits << " hits  dismissed" << endmsg;
 
@@ -326,14 +343,18 @@ std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAss
 }
 
 StatusCode DDPlanarDigi::finalize() {
-  auto file = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
+  auto file  = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
   auto names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
+  for (int i = 0; i < hSize; ++i) {
+    info() << "Standard deviation: " << m_sigma[i].get()->standard_deviation() << endmsg;
+    info() << "Mean: " << m_sigma[i].get()->mean() << endmsg;
+  }
   auto it = names.begin();
   for (auto& h : m_histograms) {
     std::string name = "";
     // Name that will appear in the stats table
     std::string    histName = *it;
-    nlohmann::json json = *h;
+    nlohmann::json json     = *h;
     auto [histo, dir] =
         Gaudi::Histograming::Sink::jsonToRootHistogram<Gaudi::Histograming::Sink::Traits<false, TH1D, 1>>(
             name, histName, json);
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 763fc38..4e0f2fd 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -19,6 +19,7 @@
 #ifndef DDPLANARDIGI_H
 #define DDPLANARDIGI_H
 
+#include <Gaudi/Accumulators.h>
 #include "Gaudi/Accumulators/Histogram.h"
 #include "Gaudi/Property.h"
 #include "GaudiAlg/Transformer.h"
@@ -124,7 +125,10 @@ struct DDPlanarDigi final
 
   const dd4hep::rec::SurfaceMap*                                        surfaceMap;
   std::array<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>, hSize> m_histograms;
-  std::string                                                           m_collName;
+  std::array<std::unique_ptr<Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>>,
+             hSize>
+              m_sigma;
+  std::string m_collName;
 
   inline static thread_local std::mt19937 m_engine;
   SmartIF<IGeoSvc>                        m_geoSvc;

From 81af3fbe675d1daeeb40b7b53a7c7720089f90d8 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Sat, 17 Feb 2024 08:44:09 +0100
Subject: [PATCH 15/24] Update to compile with the k4FWCore#173

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp |  8 ++++----
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h   | 12 ++++++------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index f37ad56..95a6cb2 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -40,11 +40,11 @@
 DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     : MultiTransformer(name, svcLoc,
                        {
-                           KeyValue("SimTrackerHitCollectionName", "SimTrackerHits"),
-                           KeyValue("HeaderName", "EventHeader"),
+                         KeyValues("SimTrackerHitCollectionName", {"SimTrackerHits"}),
+                           KeyValues("HeaderName", {"EventHeader"}),
                        },
-                       {KeyValue("TrackerHitCollectionName", "VTXTrackerHits"),
-                        KeyValue("SimTrkHitRelCollection", "VTXTrackerHitRelations")}) {
+                       {KeyValues("TrackerHitCollectionName", {"VTXTrackerHits"}),
+                        KeyValues("SimTrkHitRelCollection", {"VTXTrackerHitRelations"})}) {
   m_uidSvc = service<IUniqueIDGenSvc>("UniqueIDGenSvc", true);
   if (!m_uidSvc) {
     error() << "Unable to get UniqueIDGenSvc" << endmsg;
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 4e0f2fd..80024d9 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -24,15 +24,16 @@
 #include "Gaudi/Property.h"
 #include "GaudiAlg/Transformer.h"
 #include "GaudiKernel/RndmGenerators.h"
-#include "k4FWCore/BaseClass.h"
-#include "k4Interface/IGeoSvc.h"
-#include "k4Interface/IUniqueIDGenSvc.h"
 
 #include "edm4hep/EventHeaderCollection.h"
 #include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/TrackerHitPlaneCollection.h"
 
+#include "k4FWCore/Transformer.h"
+#include "k4Interface/IGeoSvc.h"
+#include "k4Interface/IUniqueIDGenSvc.h"
+
 #include "DDRec/SurfaceManager.h"
 
 #include <random>
@@ -74,10 +75,9 @@
 enum { hu = 0, hv, hT, hitE, hitsAccepted, diffu, diffv, diffT, hSize };
 
 struct DDPlanarDigi final
-    : Gaudi::Functional::MultiTransformer<
+    : k4FWCore::MultiTransformer<
           std::tuple<edm4hep::TrackerHitPlaneCollection, edm4hep::MCRecoTrackerHitPlaneAssociationCollection>(
-              const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&),
-          BaseClass_t> {
+              const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&)> {
   DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc);
 
   StatusCode initialize() override;

From 44ef8d67fb9b940ceead290518473aa117072c7a Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Sat, 17 Feb 2024 11:32:02 +0100
Subject: [PATCH 16/24] Update to run with k4FWCore#173

---
 .../DDPlanarDigi/components/DDPlanarDigi.cpp  | 42 +++++--------------
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h | 10 +----
 .../DDPlanarDigi/options/runDDPlanarDigi.py   | 41 +++++++++---------
 3 files changed, 32 insertions(+), 61 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 95a6cb2..1e24478 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -23,7 +23,7 @@
 #include "edm4hep/SimTrackerHit.h"
 #include "edm4hep/TrackerHitPlaneCollection.h"
 
-#include "Gaudi/Accumulators/Histogram.h"
+#include "Gaudi/Accumulators/RootHistogram.h"
 #include "Gaudi/Histograming/Sink/Utils.h"
 
 #include "DD4hep/BitFieldCoder.h"
@@ -57,27 +57,17 @@ DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     throw std::runtime_error("DDPlanarDigi: Inconsistent number of resolutions given for U and V coordinate");
   }
 
-  m_histograms[hu].reset(new Gaudi::Accumulators::Histogram<1>{this, "hu", "smearing u", {50, -5., +5.}});
-  m_histograms[hv].reset(new Gaudi::Accumulators::Histogram<1>{this, "hv", "smearing v", {50, -5., +5.}});
-  m_histograms[hT].reset(new Gaudi::Accumulators::Histogram<1>{this, "hT", "smearing time", {50, -5., +5.}});
+  m_histograms[hu].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hu", "smearing u", {50, -5., +5.}});
+  m_histograms[hv].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hv", "smearing v", {50, -5., +5.}});
+  m_histograms[hT].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hT", "smearing time", {50, -5., +5.}});
 
-  m_histograms[diffu].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffu", "diff u", {1000, -.1, +.1}});
-  m_histograms[diffv].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffv", "diff v", {1000, -.1, +.1}});
-  m_histograms[diffT].reset(new Gaudi::Accumulators::Histogram<1>{this, "diffT", "diff time", {1000, -5., +5.}});
+  m_histograms[diffu].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "diffu", "diff u", {1000, -.1, +.1}});
+  m_histograms[diffv].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "diffv", "diff v", {1000, -.1, +.1}});
+  m_histograms[diffT].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "diffT", "diff time", {1000, -5., +5.}});
 
-  m_histograms[hitE].reset(new Gaudi::Accumulators::Histogram<1>{this, "hitE", "hitEnergy in keV", {1000, 0, 200}});
+  m_histograms[hitE].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hitE", "hitEnergy in keV", {1000, 0, 200}});
   m_histograms[hitsAccepted].reset(
-      new Gaudi::Accumulators::Histogram<1>{this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
-
-  m_sigma[hu].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hv].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hT].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[diffu].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[diffv].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[diffT].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hitE].reset(new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
-  m_sigma[hitsAccepted].reset(
-      new Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>{});
+      new Gaudi::Accumulators::RootHistogram<1>{this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
 
   m_geoSvc = serviceLocator()->service(m_geoSvcName);
 }
@@ -127,7 +117,6 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
 
   for (const auto& hit : simTrackerHits) {
     ++(*m_histograms[hitE])[hit.getEDep() * (dd4hep::GeV / dd4hep::keV)];
-    *m_sigma[hitE] += hit.getEDep() * (dd4hep::GeV / dd4hep::keV);
 
     if (hit.getEDep() < m_minEnergy) {
       debug() << "Hit with insufficient energy " << hit.getEDep() * (dd4hep::GeV / dd4hep::keV) << " keV" << endmsg;
@@ -183,9 +172,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
 
       double tSmear = resT > 0 ? dist(m_engine) * resT : 0;
       ++(*m_histograms[hT])[resT > 0 ? tSmear / resT : 0];
-      *m_sigma[hT] += resT > 0 ? tSmear / resT : 0;
       ++(*m_histograms[diffT])[tSmear];
-      *m_sigma[diffT] += tSmear;
 
       hitT += tSmear;
       debug() << "smeared hit at T: " << hit.getTime() << " ns to T: " << hitT
@@ -264,14 +251,10 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
         newPos    = newPosTmp;
 
         ++(*m_histograms[hu])[uSmear / resU];
-        *m_sigma[hu] += uSmear / resU;
         ++(*m_histograms[hv])[vSmear / resV];
-        *m_sigma[hv] += vSmear / resV;
 
         ++(*m_histograms[diffu])[uSmear];
-        *m_sigma[diffu] += uSmear;
         ++(*m_histograms[diffv])[vSmear];
-        *m_sigma[diffv] += vSmear;
 
         break;
       }
@@ -335,7 +318,6 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
   // Filling the fraction of accepted hits in the event
   float accFraction = nSimHits > 0 ? float(nCreatedHits) / float(nSimHits) * 100.0 : 0.0;
   ++(*m_histograms[hitsAccepted])[accFraction];
-  *m_sigma[hitsAccepted] += accFraction;
 
   debug() << "Created " << nCreatedHits << " hits, " << nDismissedHits << " hits  dismissed" << endmsg;
 
@@ -344,11 +326,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
 
 StatusCode DDPlanarDigi::finalize() {
   auto file  = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
-  auto names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
-  for (int i = 0; i < hSize; ++i) {
-    info() << "Standard deviation: " << m_sigma[i].get()->standard_deviation() << endmsg;
-    info() << "Mean: " << m_sigma[i].get()->mean() << endmsg;
-  }
+  std::vector<const char*> names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
   auto it = names.begin();
   for (auto& h : m_histograms) {
     std::string name = "";
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 80024d9..3dca206 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -19,8 +19,7 @@
 #ifndef DDPLANARDIGI_H
 #define DDPLANARDIGI_H
 
-#include <Gaudi/Accumulators.h>
-#include "Gaudi/Accumulators/Histogram.h"
+#include "Gaudi/Accumulators/RootHistogram.h"
 #include "Gaudi/Property.h"
 #include "GaudiAlg/Transformer.h"
 #include "GaudiKernel/RndmGenerators.h"
@@ -67,8 +66,6 @@
  * (default value lcio::ILDDetID::VXD) <br>
  * <br>
  * 
- * @author F.Gaede CERN/DESY, S. Aplin DESY
- * @date Dec 2014
  * Originally in https://github.com/iLCSoft/MarlinTrkProcessors/blob/master/source/Digitisers/include/DDPlanarDigi.h
  */
 
@@ -124,10 +121,7 @@ struct DDPlanarDigi final
                                                 "Output file name for the histograms"};
 
   const dd4hep::rec::SurfaceMap*                                        surfaceMap;
-  std::array<std::unique_ptr<Gaudi::Accumulators::Histogram<1>>, hSize> m_histograms;
-  std::array<std::unique_ptr<Gaudi::Accumulators::SigmaAccumulator<Gaudi::Accumulators::atomicity::full, double>>,
-             hSize>
-              m_sigma;
+  std::array<std::unique_ptr<Gaudi::Accumulators::RootHistogram<1>>, hSize> m_histograms;
   std::string m_collName;
 
   inline static thread_local std::mt19937 m_engine;
diff --git a/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py b/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
index 16358eb..f2dc7c6 100644
--- a/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
+++ b/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
@@ -17,7 +17,8 @@
 # limitations under the License.
 #
 from Gaudi.Configuration import INFO
-from Configurables import ApplicationMgr, k4DataSvc, PodioInput, PodioOutput
+from k4FWCore import ApplicationMgr, IOSvc
+from Configurables import EventDataSvc
 from Configurables import DDPlanarDigi
 from Configurables import GeoSvc
 from Configurables import UniqueIDGenSvc
@@ -30,30 +31,28 @@
 geoservice.OutputLevel = INFO
 geoservice.EnableGeant4Geo = False
 
-processor = DDPlanarDigi()
-processor.SubDetectorName = "Vertex"
-processor.IsStrip = False
-processor.ResolutionU = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
-processor.ResolutionV = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
-processor.SimTrackerHitCollectionName = "VertexBarrelCollection"
-processor.SimTrkHitRelCollection = "VXDTrackerHitRelations"
-processor.TrackerHitCollectionName = "VXDTrackerHits"
+digi = DDPlanarDigi()
+digi.SubDetectorName = "Vertex"
+digi.IsStrip = False
+digi.ResolutionU = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
+digi.ResolutionV = [0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
+digi.SimTrackerHitCollectionName = "VertexBarrelCollection"
+digi.SimTrkHitRelCollection = "VXDTrackerHitRelations"
+digi.TrackerHitCollectionName = "VXDTrackerHits"
+digi.OutputFileName = "planar_digi_histograms.root"
 
-data_svc = k4DataSvc("EventDataSvc")
-data_svc.input = "input.root"
+iosvc = IOSvc()
+iosvc.input = "input.root"
+iosvc.output = "output_digi.root"
 
-inp = PodioInput()
-inp.collections = [
-    "VertexBarrelCollection",
-    "EventHeader",
-]
+# inp.collections = [
+#     "VertexBarrelCollection",
+#     "EventHeader",
+# ]
 
-out = PodioOutput("out")
-out.filename = "planar_digi_histograms.root"
-
-ApplicationMgr(TopAlg=[inp, processor, out],
+ApplicationMgr(TopAlg=[digi],
                EvtSel="NONE",
                EvtMax=-1,
-               ExtSvc=[data_svc],
+               ExtSvc=[EventDataSvc("EventDataSvc")],
                OutputLevel=INFO,
                )

From d5b1a099780354039614d2ebe160b73a53483172 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Thu, 29 Feb 2024 10:25:25 +0100
Subject: [PATCH 17/24] Use std::format

---
 .../DDPlanarDigi/components/DDPlanarDigi.cpp  | 21 +++++++------------
 1 file changed, 8 insertions(+), 13 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 1e24478..18fdc96 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -33,14 +33,14 @@
 #include <TFile.h>
 
 #include <cmath>
+#include <format>
 #include <iostream>
 #include <random>
-#include <sstream>
 
 DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     : MultiTransformer(name, svcLoc,
                        {
-                         KeyValues("SimTrackerHitCollectionName", {"SimTrackerHits"}),
+                           KeyValues("SimTrackerHitCollectionName", {"SimTrackerHits"}),
                            KeyValues("HeaderName", {"EventHeader"}),
                        },
                        {KeyValues("TrackerHitCollectionName", {"VTXTrackerHits"}),
@@ -66,8 +66,8 @@ DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
   m_histograms[diffT].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "diffT", "diff time", {1000, -5., +5.}});
 
   m_histograms[hitE].reset(new Gaudi::Accumulators::RootHistogram<1>{this, "hitE", "hitEnergy in keV", {1000, 0, 200}});
-  m_histograms[hitsAccepted].reset(
-      new Gaudi::Accumulators::RootHistogram<1>{this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
+  m_histograms[hitsAccepted].reset(new Gaudi::Accumulators::RootHistogram<1>{
+      this, "hitsAccepted", "Fraction of accepted hits [%]", {201, 0, 100.5}});
 
   m_geoSvc = serviceLocator()->service(m_geoSvcName);
 }
@@ -80,9 +80,7 @@ StatusCode DDPlanarDigi::initialize() {
   surfaceMap                 = surfMan->map(m_subDetName.value());
 
   if (!surfaceMap) {
-    std::stringstream err;
-    err << " Could not find surface map for detector: " << det.name() << " in SurfaceManager ";
-    throw std::runtime_error(err.str());
+    throw std::runtime_error(std::format("Could not find surface map for detector: {} in SurfaceManager", det.name()));
   }
 
   // Get and store the name for a debug message later
@@ -129,10 +127,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
     dd4hep::rec::SurfaceMap::const_iterator sI = surfaceMap->find(cellID0);
 
     if (sI == surfaceMap->end()) {
-      std::stringstream err;
-      err << " DDPlanarDigi::processEvent(): no surface found for cellID : " << std::endl;
-      // <<   cellid_decoder( hit ).valueString() ;
-      throw std::runtime_error(err.str());
+      throw std::runtime_error(std::format("DDPlanarDigi::processEvent(): no surface found for cellID : {}", cellID0));
     }
 
     const dd4hep::rec::ISurface* surf  = sI->second;
@@ -325,9 +320,9 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
 }
 
 StatusCode DDPlanarDigi::finalize() {
-  auto file  = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
+  auto                     file  = TFile::Open(m_outputFileName.value().c_str(), "RECREATE");
   std::vector<const char*> names = {"hu", "hv", "hT", "hitE", "hitsAccepted", "diffu", "diffv", "diffT", "hSize"};
-  auto it = names.begin();
+  auto                     it    = names.begin();
   for (auto& h : m_histograms) {
     std::string name = "";
     // Name that will appear in the stats table

From 89eec48950f293ad1db0f6116e2e46cb551987d9 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Thu, 29 Feb 2024 12:45:38 +0100
Subject: [PATCH 18/24] Use fmtlib

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 18fdc96..184f078 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -33,9 +33,9 @@
 #include <TFile.h>
 
 #include <cmath>
-#include <format>
 #include <iostream>
 #include <random>
+#include <fmt/format.h>
 
 DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     : MultiTransformer(name, svcLoc,
@@ -80,7 +80,7 @@ StatusCode DDPlanarDigi::initialize() {
   surfaceMap                 = surfMan->map(m_subDetName.value());
 
   if (!surfaceMap) {
-    throw std::runtime_error(std::format("Could not find surface map for detector: {} in SurfaceManager", det.name()));
+    throw std::runtime_error(fmt::format("Could not find surface map for detector: {} in SurfaceManager", det.name()));
   }
 
   // Get and store the name for a debug message later
@@ -127,7 +127,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
     dd4hep::rec::SurfaceMap::const_iterator sI = surfaceMap->find(cellID0);
 
     if (sI == surfaceMap->end()) {
-      throw std::runtime_error(std::format("DDPlanarDigi::processEvent(): no surface found for cellID : {}", cellID0));
+      throw std::runtime_error(fmt::format("DDPlanarDigi::processEvent(): no surface found for cellID : {}", cellID0));
     }
 
     const dd4hep::rec::ISurface* surf  = sI->second;

From 867e269521c709926106e381bbb1e062b154f670 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Mon, 24 Jun 2024 17:58:29 +0200
Subject: [PATCH 19/24] Remove a few includes

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 1 -
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h   | 2 --
 2 files changed, 3 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 184f078..79df11a 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -33,7 +33,6 @@
 #include <TFile.h>
 
 #include <cmath>
-#include <iostream>
 #include <random>
 #include <fmt/format.h>
 
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 3dca206..c9c0df8 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -21,8 +21,6 @@
 
 #include "Gaudi/Accumulators/RootHistogram.h"
 #include "Gaudi/Property.h"
-#include "GaudiAlg/Transformer.h"
-#include "GaudiKernel/RndmGenerators.h"
 
 #include "edm4hep/EventHeaderCollection.h"
 #include "edm4hep/MCRecoTrackerHitPlaneAssociationCollection.h"

From c270a61d4cca0e2de55d1679808be9ae80cfcf6e Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Mon, 24 Jun 2024 18:06:34 +0200
Subject: [PATCH 20/24] Fix pre-commit

---
 .../DDPlanarDigi/components/DDPlanarDigi.cpp   |  4 ++--
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h  | 18 +++++++++---------
 k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py |  2 +-
 test/CMakeLists.txt                            |  2 +-
 4 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 79df11a..127a8ed 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2014-2023 Key4hep-Project.
+ * Copyright (c) 2020-2024 Key4hep-Project.
  *
  * This file is part of Key4hep.
  * See https://key4hep.github.io/key4hep-doc/ for further info.
@@ -32,9 +32,9 @@
 
 #include <TFile.h>
 
+#include <fmt/format.h>
 #include <cmath>
 #include <random>
-#include <fmt/format.h>
 
 DDPlanarDigi::DDPlanarDigi(const std::string& name, ISvcLocator* svcLoc)
     : MultiTransformer(name, svcLoc,
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index c9c0df8..8a48ac3 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2014-2023 Key4hep-Project.
+ * Copyright (c) 2020-2024 Key4hep-Project.
  *
  * This file is part of Key4hep.
  * See https://key4hep.github.io/key4hep-doc/ for further info.
@@ -38,13 +38,13 @@
 #include <vector>
 
 /** ======= DDPlanarDigi ========== <br>
- * Creates TrackerHits from SimTrackerHits, smearing them according to the input parameters. 
+ * Creates TrackerHits from SimTrackerHits, smearing them according to the input parameters.
  * The positions of "digitized" TrackerHits are obtained by gaussian smearing positions
- * of SimTrackerHits perpendicular and along the ladder according to the specified point resolutions. 
+ * of SimTrackerHits perpendicular and along the ladder according to the specified point resolutions.
  * The geometry of the surface is retreived from DDRec::Surface associated to the hit via cellID.
- * 
- * 
- * <h4>Input collections and prerequisites</h4> 
+ *
+ *
+ * <h4>Input collections and prerequisites</h4>
  * Processor requires a collection of SimTrackerHits <br>
  * <h4>Output</h4>
  * Processor produces collection of smeared TrackerHits<br>
@@ -63,7 +63,7 @@
  * @param Sub_Detector_ID ID of Sub-Detector using UTIL/ILDConf.h from lcio <br>
  * (default value lcio::ILDDetID::VXD) <br>
  * <br>
- * 
+ *
  * Originally in https://github.com/iLCSoft/MarlinTrkProcessors/blob/master/source/Digitisers/include/DDPlanarDigi.h
  */
 
@@ -118,9 +118,9 @@ struct DDPlanarDigi final
   Gaudi::Property<std::string> m_outputFileName{this, "OutputFileName", "planar_digi_histograms.root",
                                                 "Output file name for the histograms"};
 
-  const dd4hep::rec::SurfaceMap*                                        surfaceMap;
+  const dd4hep::rec::SurfaceMap*                                            surfaceMap;
   std::array<std::unique_ptr<Gaudi::Accumulators::RootHistogram<1>>, hSize> m_histograms;
-  std::string m_collName;
+  std::string                                                               m_collName;
 
   inline static thread_local std::mt19937 m_engine;
   SmartIF<IGeoSvc>                        m_geoSvc;
diff --git a/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py b/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
index f2dc7c6..75d58c6 100644
--- a/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
+++ b/k4Reco/DDPlanarDigi/options/runDDPlanarDigi.py
@@ -1,5 +1,5 @@
 #
-# Copyright (c) 2020-2023 Key4hep-Project.
+# Copyright (c) 2020-2024 Key4hep-Project.
 #
 # This file is part of Key4hep.
 # See https://key4hep.github.io/key4hep-doc/ for further info.
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index a2d02c8..3527b6b 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,5 +1,5 @@
 #[[
-Copyright (c) 2019-2024 Key4hep-Project.
+Copyright (c) 2020-2024 Key4hep-Project.
 
 This file is part of Key4hep.
 See https://key4hep.github.io/key4hep-doc/ for further info.

From b90afbcb2230180d6435eeed25fccf4debf3611e Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Mon, 24 Jun 2024 18:39:05 +0200
Subject: [PATCH 21/24] Use TRandom2

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 11 ++++++-----
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h   |  4 +++-
 2 files changed, 9 insertions(+), 6 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index 127a8ed..c371954 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -30,7 +30,8 @@
 #include "DD4hep/DD4hepUnits.h"
 #include "DD4hep/Detector.h"
 
-#include <TFile.h>
+#include "TFile.h"
+#include "TMath.h"
 
 #include <fmt/format.h>
 #include <cmath>
@@ -94,7 +95,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
   auto seed = m_uidSvc->getUniqueID(headers[0].getEventNumber(), headers[0].getRunNumber(), this->name());
   debug() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
           << headers[0].getRunNumber() << endmsg;
-  m_engine.seed(seed);
+  m_engine.SetSeed(seed);
   // For rng calls, use the fact that drawing from a
   // Gaussian with mean mu and variance sigma^2 is the same as drawing from
   // a normal distribution and then multiplying by sigma and adding mu
@@ -164,7 +165,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
     if (m_resTLayer.size() and m_resTLayer[0] > 0) {
       float resT = m_resTLayer.size() > 1 ? m_resTLayer[layer] : m_resTLayer[0];
 
-      double tSmear = resT > 0 ? dist(m_engine) * resT : 0;
+      double tSmear = resT > 0 ? m_engine.Gaus(0, 1) * resT : 0;
       ++(*m_histograms[hT])[resT > 0 ? tSmear / resT : 0];
       ++(*m_histograms[diffT])[tSmear];
 
@@ -213,8 +214,8 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
     while (tries < m_maxTries) {
       // if( tries > 0 ) debug() << "retry smearing for " <<  cellid_decoder( hit ).valueString() << " : retries " << tries << endmsg;
 
-      double uSmear = dist(m_engine) * resU;
-      double vSmear = dist(m_engine) * resV;
+      double uSmear = m_engine.Gaus(0, 1) * resU;
+      double vSmear = m_engine.Gaus(0, 1) * resV;
 
       dd4hep::rec::Vector3D newPosTmp;
       if (m_isStrip) {
diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 8a48ac3..6f3a55d 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -33,6 +33,8 @@
 
 #include "DDRec/SurfaceManager.h"
 
+#include "TRandom2.h"
+
 #include <random>
 #include <string>
 #include <vector>
@@ -122,7 +124,7 @@ struct DDPlanarDigi final
   std::array<std::unique_ptr<Gaudi::Accumulators::RootHistogram<1>>, hSize> m_histograms;
   std::string                                                               m_collName;
 
-  inline static thread_local std::mt19937 m_engine;
+  inline static thread_local TRandom2     m_engine;
   SmartIF<IGeoSvc>                        m_geoSvc;
   SmartIF<IUniqueIDGenSvc>                m_uidSvc;
 };

From 01cbb771e84e37607df8c55674232b03e765ef6b Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Mon, 24 Jun 2024 18:39:21 +0200
Subject: [PATCH 22/24] Fix pre-commit

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.h | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
index 6f3a55d..4e5f0ce 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.h
@@ -124,9 +124,9 @@ struct DDPlanarDigi final
   std::array<std::unique_ptr<Gaudi::Accumulators::RootHistogram<1>>, hSize> m_histograms;
   std::string                                                               m_collName;
 
-  inline static thread_local TRandom2     m_engine;
-  SmartIF<IGeoSvc>                        m_geoSvc;
-  SmartIF<IUniqueIDGenSvc>                m_uidSvc;
+  inline static thread_local TRandom2 m_engine;
+  SmartIF<IGeoSvc>                    m_geoSvc;
+  SmartIF<IUniqueIDGenSvc>            m_uidSvc;
 };
 
 DECLARE_COMPONENT(DDPlanarDigi)

From beb2f55933a524f8b8d5f89f59dbd8278db4dd2d Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Mon, 24 Jun 2024 18:51:29 +0200
Subject: [PATCH 23/24] Link to ROOT::Core and ROOT::MathCore

---
 k4Reco/CMakeLists.txt | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/k4Reco/CMakeLists.txt b/k4Reco/CMakeLists.txt
index 3c35720..14e7d25 100644
--- a/k4Reco/CMakeLists.txt
+++ b/k4Reco/CMakeLists.txt
@@ -28,6 +28,8 @@ gaudi_add_module(k4RecoPlugins
                       EDM4HEP::edm4hep
                       DD4hep::DDRec
                       DD4hep::DDCore
+                      ROOT::Core
+                      ROOT::MathCore
                       )
 target_include_directories(k4RecoPlugins PUBLIC
   $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}/include>

From c01b8f7af25ce4bcdb7c6214eaad12acb0e72550 Mon Sep 17 00:00:00 2001
From: jmcarcell <jmcarcell@users.noreply.github.com>
Date: Mon, 24 Jun 2024 19:02:14 +0200
Subject: [PATCH 24/24] Remove a comment and use the Gaus function with the
 right sigma

---
 k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
index c371954..b7dfce9 100644
--- a/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
+++ b/k4Reco/DDPlanarDigi/components/DDPlanarDigi.cpp
@@ -96,10 +96,6 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
   debug() << "Using seed " << seed << " for event " << headers[0].getEventNumber() << " and run "
           << headers[0].getRunNumber() << endmsg;
   m_engine.SetSeed(seed);
-  // For rng calls, use the fact that drawing from a
-  // Gaussian with mean mu and variance sigma^2 is the same as drawing from
-  // a normal distribution and then multiplying by sigma and adding mu
-  auto dist = std::normal_distribution<double>(0, 1);
 
   int nCreatedHits   = 0;
   int nDismissedHits = 0;
@@ -165,7 +161,7 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
     if (m_resTLayer.size() and m_resTLayer[0] > 0) {
       float resT = m_resTLayer.size() > 1 ? m_resTLayer[layer] : m_resTLayer[0];
 
-      double tSmear = resT > 0 ? m_engine.Gaus(0, 1) * resT : 0;
+      double tSmear = resT > 0 ? m_engine.Gaus(0, resT) : 0;
       ++(*m_histograms[hT])[resT > 0 ? tSmear / resT : 0];
       ++(*m_histograms[diffT])[tSmear];
 
@@ -214,8 +210,8 @@ DDPlanarDigi::operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits,
     while (tries < m_maxTries) {
       // if( tries > 0 ) debug() << "retry smearing for " <<  cellid_decoder( hit ).valueString() << " : retries " << tries << endmsg;
 
-      double uSmear = m_engine.Gaus(0, 1) * resU;
-      double vSmear = m_engine.Gaus(0, 1) * resV;
+      double uSmear = m_engine.Gaus(0, resU);
+      double vSmear = m_engine.Gaus(0, resV);
 
       dd4hep::rec::Vector3D newPosTmp;
       if (m_isStrip) {