From a999d24ad1a9a65f82c93a84a11d1db553f932a8 Mon Sep 17 00:00:00 2001 From: Brieuc Francois Date: Thu, 8 Aug 2024 09:42:11 +0200 Subject: [PATCH] Transformer to create tracks from gen particles (#21) * Transformer to create tracks from gen particles * Move away from BaseClass_t * Adapt to the new way Transformer work * Migrate to IOSvc * [PlotTrackHitResiduals] Remove overflow bin and use RootHistSvc * [TracksFromGenParticles] Fix tests * [PlotTrackHitResiduals] Set histogram name * Fix tests following removal of GaudiAlg * Fix tests following removal of GaudiAlg --- ARCdigi/test/runARCdigitizer.py | 5 -- DCHdigi/test/runDCHsimpleDigitizer.py | 5 -- .../test/runDCHsimpleDigitizerExtendedEdm.py | 5 -- Tracking/CMakeLists.txt | 8 ++ Tracking/components/PlotTrackHitResiduals.cpp | 81 +++++++++++++++++ .../components/TracksFromGenParticles.cpp | 90 +++++++++++++++++++ ...nGenFitTrackingOnSimplifiedDriftChamber.py | 5 -- Tracking/test/runTracksFromGenParticles.py | 50 +++++++++++ VTXdigi/test/runVTXdigitizer.py | 6 -- 9 files changed, 229 insertions(+), 26 deletions(-) create mode 100644 Tracking/components/PlotTrackHitResiduals.cpp create mode 100644 Tracking/components/TracksFromGenParticles.cpp create mode 100644 Tracking/test/runTracksFromGenParticles.py diff --git a/ARCdigi/test/runARCdigitizer.py b/ARCdigi/test/runARCdigitizer.py index 712a3bb..9ecc734 100644 --- a/ARCdigi/test/runARCdigitizer.py +++ b/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 diff --git a/DCHdigi/test/runDCHsimpleDigitizer.py b/DCHdigi/test/runDCHsimpleDigitizer.py index e16b382..d310c90 100644 --- a/DCHdigi/test/runDCHsimpleDigitizer.py +++ b/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, diff --git a/DCHdigi/test/runDCHsimpleDigitizerExtendedEdm.py b/DCHdigi/test/runDCHsimpleDigitizerExtendedEdm.py index 3d30271..4990c46 100644 --- a/DCHdigi/test/runDCHsimpleDigitizerExtendedEdm.py +++ b/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, diff --git a/Tracking/CMakeLists.txt b/Tracking/CMakeLists.txt index f441d5c..617d1eb 100644 --- a/Tracking/CMakeLists.txt +++ b/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 $ $ + ${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() diff --git a/Tracking/components/PlotTrackHitResiduals.cpp b/Tracking/components/PlotTrackHitResiduals.cpp new file mode 100644 index 0000000..a32e6b9 --- /dev/null +++ b/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 + +// ROOT +#include "TH1D.h" + +// k4FWCore +#include "k4FWCore/Consumer.h" + +#include + +/** @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 { + 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 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) diff --git a/Tracking/components/TracksFromGenParticles.cpp b/Tracking/components/TracksFromGenParticles.cpp new file mode 100644 index 0000000..af03cb8 --- /dev/null +++ b/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 + +// k4FWCore +#include "k4FWCore/Transformer.h" + +#include + +/** @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(const edm4hep::MCParticleCollection&)> { + TracksFromGenParticles(const std::string& name, ISvcLocator* svcLoc) + : MultiTransformer( + name, svcLoc, + {KeyValues("InputGenParticles", {"MCParticles"})}, + {KeyValues("OutputTracks", {"TracksFromGenParticles"}), + KeyValues("OutputMCRecoTrackParticleAssociation", {"TracksFromGenParticlesAssociation"})}) { + } + +std::tuple 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 m_Bz{this, "Bz", 2., "Z component of the (assumed constant) magnetic field in Tesla."}; +}; + +DECLARE_COMPONENT(TracksFromGenParticles) diff --git a/Tracking/test/runGenFitTrackingOnSimplifiedDriftChamber.py b/Tracking/test/runGenFitTrackingOnSimplifiedDriftChamber.py index e726835..1427308 100644 --- a/Tracking/test/runGenFitTrackingOnSimplifiedDriftChamber.py +++ b/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, diff --git a/Tracking/test/runTracksFromGenParticles.py b/Tracking/test/runTracksFromGenParticles.py new file mode 100644 index 0000000..8044b80 --- /dev/null +++ b/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, +) diff --git a/VTXdigi/test/runVTXdigitizer.py b/VTXdigi/test/runVTXdigitizer.py index b217704..49d2031 100644 --- a/VTXdigi/test/runVTXdigitizer.py +++ b/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,