Transformer to create tracks from gen particles

ARCdigi/test/runARCdigitizer.py

@@ -34,13 +34,8 @@
 arc_digitizer.AuditExecute = True
 podiooutput.AuditExecute = True
 
-from Configurables import EventCounter
-event_counter = EventCounter('event_counter')
-event_counter.Frequency = 1
-
 ApplicationMgr(
     TopAlg = [
-        event_counter,
         podioinput,
         arc_digitizer,
         podiooutput

DCHdigi/test/runDCHsimpleDigitizer.py

@@ -129,14 +129,9 @@
 geantsim.AuditExecute = True
 out.AuditExecute = True
 
-from Configurables import EventCounter
-event_counter = EventCounter('event_counter')
-event_counter.Frequency = 1
-
 from Configurables import ApplicationMgr
 ApplicationMgr(
     TopAlg = [
-              event_counter,
               genAlg,
               hepmc_converter,
               geantsim,

DCHdigi/test/runDCHsimpleDigitizerExtendedEdm.py

@@ -162,14 +162,9 @@
 geantsim.AuditExecute = True
 out.AuditExecute = True
 
-from Configurables import EventCounter
-event_counter = EventCounter('event_counter')
-event_counter.Frequency = 1
-
 from Configurables import ApplicationMgr
 ApplicationMgr(
     TopAlg = [
-              event_counter,
               genAlg,
               hepmc_converter,
               geantsim,

Tracking/CMakeLists.txt

@@ -3,11 +3,13 @@ set(PackageName Tracking)
 project(${PackageName})
 
 #find_package(GenFit)
+FIND_PACKAGE(MarlinUtil)
 
 #if (GenFit_FOUND)
 
 file(GLOB sources
     ${PROJECT_SOURCE_DIR}/src/*.cpp
+    ${PROJECT_SOURCE_DIR}/components/*.cpp
 )
 
 file(GLOB headers
@@ -26,6 +28,7 @@ gaudi_add_module(${PackageName}
 target_include_directories(${PackageName} PUBLIC
   $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
   $<INSTALL_INTERFACE:include>
+  ${MarlinUtil_INCLUDE_DIRS}
 )
 
 set_target_properties(${PackageName} PROPERTIES PUBLIC_HEADER "${headers}")
@@ -44,4 +47,9 @@ install(TARGETS ${PackageName}
 )
 
 install(FILES ${scripts} DESTINATION test)
+
+SET(test_name "test_TracksFromGenParticles")
+ADD_TEST(NAME ${test_name} COMMAND k4run test/runTracksFromGenParticles.py)
+set_test_env(${test_name})
+
 #endif()

Tracking/components/PlotTrackHitResiduals.cpp

@@ -0,0 +1,81 @@
+// Gaudi
+#include "Gaudi/Property.h"
+#include "Gaudi/Accumulators/RootHistogram.h"
+#include "Gaudi/Histograming/Sink/Utils.h"
+
+// edm4hep
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+
+// marlin
+#include <marlinutil/HelixClass_double.h>
+
+// ROOT
+#include "TH1D.h"
+
+// k4FWCore
+#include "k4FWCore/Consumer.h"
+
+#include <string>
+
+/** @class PlotTrackHitDistances
+ *
+ *  Gaudi consumer that generates a residual distribution (mm) by comparing the helix from Track AtIP and simHit position.
+ *  This is intended to be used on tracks produced from gen particles i.e. which do not have real hits attached to them.
+ *
+ *  @author Brieuc Francois
+ */
+
+struct PlotTrackHitDistances final
+  : k4FWCore::Consumer<void(const edm4hep::SimTrackerHitCollection&, const edm4hep::MCRecoTrackParticleAssociationCollection&)> {
+  PlotTrackHitDistances(const std::string& name, ISvcLocator* svcLoc)
+      : Consumer(
+            name, svcLoc,
+            {
+            KeyValues("InputSimTrackerHits", {"DCHCollection"}),
+            KeyValues("InputTracksFromGenParticlesAssociation", {"TracksFromGenParticlesAssociation"}),
+            }) {}
+
+  void operator()(const edm4hep::SimTrackerHitCollection& simTrackerHits, const edm4hep::MCRecoTrackParticleAssociationCollection& trackParticleAssociations) const override {
+
+    for (const auto& trackParticleAssociation : trackParticleAssociations) {
+      auto genParticle = trackParticleAssociation.getSim();
+      auto track = trackParticleAssociation.getRec();
+      edm4hep::TrackState trackStateAtIP;
+      bool found_trackStateAtIP = false;
+      for (const auto& trackState : track.getTrackStates()) {
+        if (trackState.location == edm4hep::TrackState::AtIP) {
+          trackStateAtIP = trackState;
+          found_trackStateAtIP = true;
+          break;
+        }
+      }
+      if (!found_trackStateAtIP)
+        throw std::runtime_error("No track state defined AtIP, exiting!");
+
+      // Build an helix out of the trackState
+      auto helixFromTrack = HelixClass_double();
+      helixFromTrack.Initialize_Canonical(trackStateAtIP.phi, trackStateAtIP.D0, trackStateAtIP.Z0, trackStateAtIP.omega, trackStateAtIP.tanLambda, m_Bz);
+
+      // Fill the histogram with residuals for hits attached to the same gen particle
+      for (const auto& simTrackerHit : simTrackerHits) {
+        auto simTrackerHitgenParticle = simTrackerHit.getParticle();
+        if (simTrackerHitgenParticle.getObjectID() == genParticle.getObjectID()) {
+          double simTrackerHitPosition[] = {simTrackerHit.x(), simTrackerHit.y(), simTrackerHit.z()};
+          double distances[3];
+          helixFromTrack.getDistanceToPoint(simTrackerHitPosition, distances);
+          // Distance[0] - distance in R-Phi plane, Distance[1] - distance along Z axis, Distance[2] - 3D distance
+          ++m_residualHist[distances[2]];
+        }
+      }
+    }
+    return;
+  }
+  Gaudi::Property<float> m_Bz{this, "Bz", 2., "Z component of the (assumed constant) magnetic field in Tesla."};
+  mutable Gaudi::Accumulators::Histogram<1> m_residualHist{this, "track_hits_distance_closest_approach", "Track-hit Distances", {100, 0, 1, "Distance [mm];Entries"}};
+
+};
+
+DECLARE_COMPONENT(PlotTrackHitDistances)

Tracking/components/TracksFromGenParticles.cpp

@@ -0,0 +1,90 @@
+#include "Gaudi/Property.h"
+
+// edm4hep
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+
+// marlin
+#include <marlinutil/HelixClass_double.h>
+
+// k4FWCore
+#include "k4FWCore/Transformer.h"
+
+#include <string>
+
+/** @class TracksFromGenParticles
+ *
+ *  Gaudi transformer that builds an edm4hep::TrackCollection out of an edm4hep::MCParticleCollection.
+ *  It just builds an helix out of the genParticle position, momentum, charge and user defined z component of the (constant) magnetic field.
+ *  From this helix, different edm4hep::TrackStates (AtIP, AtFirstHit, AtLastHit and AtCalorimeter) are defined. #FIXME for now these trackstates are dummy (copies of the same helix parameters)
+ *  This is meant to enable technical development needing edm4hep::Track and performance studies where having generator based trackis is a reasonable approximation.
+ *  Possible inprovement:
+ *    - Retrieve magnetic field from geometry: const DD4hep::Field::MagneticField& magneticField = detector.field(); DD4hep::DDRec::Vector3D field = magneticField.magneticField(point);
+ *    - Properly define different trackStates
+ *
+ *  @author Brieuc Francois
+ */
+
+struct TracksFromGenParticles final
+  : k4FWCore::MultiTransformer<std::tuple<edm4hep::TrackCollection, edm4hep::MCRecoTrackParticleAssociationCollection>(const edm4hep::MCParticleCollection&)> {
+  TracksFromGenParticles(const std::string& name, ISvcLocator* svcLoc)
+      : MultiTransformer(
+            name, svcLoc,
+            {KeyValues("InputGenParticles", {"MCParticles"})},
+            {KeyValues("OutputTracks", {"TracksFromGenParticles"}),
+            KeyValues("OutputMCRecoTrackParticleAssociation", {"TracksFromGenParticlesAssociation"})}) {
+  }
+
+std::tuple<edm4hep::TrackCollection, edm4hep::MCRecoTrackParticleAssociationCollection> operator()(const edm4hep::MCParticleCollection& genParticleColl) const override {
+
+    auto outputTrackCollection = edm4hep::TrackCollection();
+    auto MCRecoTrackParticleAssociationCollection = edm4hep::MCRecoTrackParticleAssociationCollection();
+
+    for (const auto& genParticle : genParticleColl) {
+      debug() << "Particle decayed in tracker: " << genParticle.isDecayedInTracker() << endmsg;
+      debug() << genParticle << endmsg;
+
+      // Building an helix out of MCParticle properties and B field
+      auto helixFromGenParticle = HelixClass_double();
+      double genParticleVertex[] = {genParticle.getVertex().x, genParticle.getVertex().y, genParticle.getVertex().z};
+      double genParticleMomentum[] = {genParticle.getMomentum().x, genParticle.getMomentum().y, genParticle.getMomentum().z};
+      helixFromGenParticle.Initialize_VP(genParticleVertex, genParticleMomentum, genParticle.getCharge(), m_Bz);
+
+      // Setting the track and trackStates properties
+      // #FIXME for now, the different trackStates are dummy
+      auto trackFromGen = edm4hep::MutableTrack();
+      auto trackState_IP = edm4hep::TrackState {};
+      trackState_IP.location = edm4hep::TrackState::AtIP;
+      trackState_IP.D0 = helixFromGenParticle.getD0();
+      trackState_IP.phi = helixFromGenParticle.getPhi0();
+      trackState_IP.omega = helixFromGenParticle.getOmega();
+      trackState_IP.Z0 = helixFromGenParticle.getZ0();
+      trackState_IP.tanLambda = helixFromGenParticle.getTanLambda();
+      trackFromGen.addToTrackStates(trackState_IP);
+      auto trackState_AtFirstHit = edm4hep::TrackState(trackState_IP);
+      trackState_AtFirstHit.location = edm4hep::TrackState::AtFirstHit;
+      trackFromGen.addToTrackStates(trackState_AtFirstHit);
+      auto trackState_AtLastHit = edm4hep::TrackState(trackState_IP);
+      trackState_AtFirstHit.location = edm4hep::TrackState::AtLastHit;
+      trackFromGen.addToTrackStates(trackState_AtLastHit);
+      auto trackState_AtCalorimeter = edm4hep::TrackState(trackState_IP);
+      trackState_AtFirstHit.location = edm4hep::TrackState::AtCalorimeter;
+      trackFromGen.addToTrackStates(trackState_AtCalorimeter);
+
+      //debug() << trackFromGen << endmsg;
+      outputTrackCollection.push_back(trackFromGen);
+
+      // Building the association between tracks and genParticles
+      auto MCRecoTrackParticleAssociation = edm4hep::MutableMCRecoTrackParticleAssociation();
+      MCRecoTrackParticleAssociation.setRec(trackFromGen);
+      MCRecoTrackParticleAssociation.setSim(genParticle);
+      MCRecoTrackParticleAssociationCollection.push_back(MCRecoTrackParticleAssociation);
+    }
+    return std::make_tuple(std::move(outputTrackCollection), std::move(MCRecoTrackParticleAssociationCollection));
+  }
+
+  Gaudi::Property<float> m_Bz{this, "Bz", 2., "Z component of the (assumed constant) magnetic field in Tesla."};
+};
+
+DECLARE_COMPONENT(TracksFromGenParticles)

Tracking/test/runGenFitTrackingOnSimplifiedDriftChamber.py

@@ -163,14 +163,9 @@
 geantsim.AuditExecute = True
 out.AuditExecute = True
 
-from Configurables import EventCounter
-event_counter = EventCounter('event_counter')
-event_counter.Frequency = 1
-
 from Configurables import ApplicationMgr
 ApplicationMgr(
     TopAlg = [
-              event_counter,
               genAlg,
               hepmc_converter,
               geantsim,

Tracking/test/runTracksFromGenParticles.py

@@ -0,0 +1,50 @@
+from k4FWCore import ApplicationMgr
+from Gaudi.Configuration import INFO, WARNING, DEBUG
+import os
+
+if not os.path.isfile("ddsim_output_edm4hep.root"):
+    os.system("ddsim --enableGun --gun.distribution uniform --gun.energy '10*GeV' --gun.particle e- --numberOfEvents 100 --outputFile ddsim_output_edm4hep.root --random.enableEventSeed --random.seed 42 --compactFile $K4GEO/FCCee/IDEA/compact/IDEA_o1_v02/IDEA_o1_v02.xml")
+
+# Loading the output of the SIM step
+from k4FWCore import IOSvc
+io_svc = IOSvc("IOSvc")
+io_svc.input = "ddsim_output_edm4hep.root"
+io_svc.output = "tracks_from_genParticle_output.root"
+
+# Calling TracksFromGenParticles
+from Configurables import TracksFromGenParticles
+tracksFromGenParticles = TracksFromGenParticles("TracksFromGenParticles",
+                                               InputGenParticles = ["MCParticles"],
+                                               OutputTracks = ["TracksFromGenParticles"],
+                                               OutputMCRecoTrackParticleAssociation = ["TracksFromGenParticlesAssociation"],
+                                               Bz = 2.0,
+                                               OutputLevel = INFO)
+
+# produce a TH1 with distances between tracks and simTrackerHits
+from Configurables import PlotTrackHitDistances, RootHistSvc
+from Configurables import Gaudi__Histograming__Sink__Root as RootHistoSink
+plotTrackHitDistances = PlotTrackHitDistances("PlotTrackHitDistances",
+                                             InputSimTrackerHits = ["CDCHHits"],
+                                             InputTracksFromGenParticlesAssociation = tracksFromGenParticles.OutputMCRecoTrackParticleAssociation, 
+                                             Bz = 2.0)
+
+hps = RootHistSvc("HistogramPersistencySvc")
+root_hist_svc = RootHistoSink("RootHistoSink")
+root_hist_svc.FileName = "TrackHitDistances.root"
+
+# Set auditor service
+from Configurables import AuditorSvc, ChronoAuditor
+chra = ChronoAuditor()
+audsvc = AuditorSvc()
+audsvc.Auditors = [chra]
+tracksFromGenParticles.AuditExecute = True
+plotTrackHitDistances.AuditExecute = True
+
+from Configurables import EventDataSvc
+ApplicationMgr(
+    TopAlg= [tracksFromGenParticles, plotTrackHitDistances],
+    EvtSel='NONE',
+    EvtMax=-1,
+    ExtSvc=[root_hist_svc, EventDataSvc("EventDataSvc"), audsvc],
+    StopOnSignal=True,
+)

VTXdigi/test/runVTXdigitizer.py

@@ -268,16 +268,11 @@
 geantsim.AuditExecute = True
 out.AuditExecute = True
 
-from Configurables import EventCounter
-event_counter = EventCounter('event_counter')
-event_counter.Frequency = 1
-
 from Configurables import ApplicationMgr
 
 # # CLD
 # ApplicationMgr(
 #     TopAlg = [
-#               event_counter,
 #               genAlg,
 #               hepmc_converter,
 #               geantsim,
@@ -293,7 +288,6 @@
 # IDEA
 ApplicationMgr(
     TopAlg = [
-              event_counter,
               genAlg,
               hepmc_converter,
               geantsim,