generated from key4hep/k4-project-template
-
Notifications
You must be signed in to change notification settings - Fork 13
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
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
- Loading branch information
Showing
9 changed files
with
229 additions
and
26 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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, | ||
| ) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters