Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Changes to run pandora on tracks created from generator-level particles #43

Merged
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 14 additions & 2 deletions Tracking/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -4,7 +4,15 @@ project(${PackageName})

#find_package(GenFit)
find_package(MarlinUtil REQUIRED)

list(APPEND CMAKE_MODULE_PATH $ENV{PANDORAPFA}/cmakemodules)
find_package(PandoraSDK REQUIRED)
FOREACH( pkg MarlinUtil PandoraSDK )
IF( ${pkg}_FOUND )
INCLUDE_DIRECTORIES( ${${pkg}_INCLUDE_DIRS} )
LINK_LIBRARIES( ${${pkg}_LIBRARIES} )
ADD_DEFINITIONS ( ${${pkg}_DEFINITIONS} )
ENDIF()
ENDFOREACH()
#if (GenFit_FOUND)

file(GLOB sources
@@ -43,7 +51,7 @@ install(TARGETS ${PackageName}
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
PUBLIC_HEADER DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}/${CMAKE_PROJECT_NAME}" COMPONENT dev
)

install(FILES ${scripts} DESTINATION test)
@@ -52,4 +60,8 @@ SET(test_name "test_TracksFromGenParticles")
ADD_TEST(NAME ${test_name} COMMAND k4run test/runTracksFromGenParticles.py)
set_test_env(${test_name})

SET(test_name "test_TracksFromGenParticlesAlg")
ADD_TEST(NAME ${test_name} COMMAND k4run test/runTracksFromGenParticlesAlg.py)
set_test_env(${test_name})

#endif()
376 changes: 376 additions & 0 deletions Tracking/components/TracksFromGenParticlesWithECalExtrapAlg.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,376 @@
// Gaudi
#include "Gaudi/Property.h"

// edm4hep
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/TrackCollection.h"
#include "edm4hep/TrackMCParticleLinkCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"

// marlin
#include <marlinutil/HelixClass_double.h>

// pandora
#include <Objects/Helix.h>

// k4FWCore
#include "k4FWCore/DataHandle.h"

// Gaudi
#include "Gaudi/Algorithm.h"
#include "GaudiKernel/ToolHandle.h"

// DD4HEP
#include <DDRec/DetectorData.h>
#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/DetType.h"
#include "DD4hep/DetectorSelector.h"

// C++
#include <string>

/** @class TracksFromGenParticlesWithECalExtrapAlg
*
* GaudiAlg version of TracksFromGenParticles, that builds an edm4hep::TrackCollection out of an edm4hep::MCParticleCollection.
* It just builds an helix out of the genParticle position, momentum, charge and z component of the (constant) magnetic field, retrieved from the detector.
* From this helix, different edm4hep::TrackStates (AtIP, AtFirstHit, AtLastHit) are defined.
* The first and last hits are defined as those with smallest and largest time in the input SimTrackerHit collections
* The algorithm also performs extrapolation to the EM calorimeter inner face is done using the positions of the barrel and endcap retrieved from the detector data extensions.
* This is meant to enable technical development needing edm4hep::Track and performance studies where having generator based tracks is a reasonable approximation.
* @author Brieuc Francois
* @author Archil Durglishvili
* @author Giovanni Marchiori
*/

class TracksFromGenParticlesWithECalExtrapAlg : public Gaudi::Algorithm {

public:
TracksFromGenParticlesWithECalExtrapAlg(const std::string& name, ISvcLocator* svcLoc);
StatusCode initialize();
StatusCode execute(const EventContext&) const;
StatusCode finalize();

private:
/// Handle for input MC particles
mutable DataHandle<edm4hep::MCParticleCollection> m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this};
/// List of input sim tracker hit collections
Gaudi::Property<std::vector<std::string>> m_inputSimTrackerHitCollections{this, "InputSimTrackerHits", {}, "Names of SimTrackerHit collections to read"};
/// the vector of input DataHandles for the tracker hit collections
std::vector<DataObjectHandleBase*> m_inputSimTrackerHitCollectionHandles;
/// Solenoid magnetic field, to be retrieved from detector
float m_Bz;
/// ECAL barrel and endcap extent
float m_eCalBarrelInnerR;
float m_eCalBarrelMaxZ;
float m_eCalEndCapInnerR;
float m_eCalEndCapOuterR;
float m_eCalEndCapInnerZ;
float m_eCalEndCapOuterZ;
/// Handle for the output track collection
mutable DataHandle<edm4hep::TrackCollection> m_tracks{"TracksFromGenParticlesAlg", Gaudi::DataHandle::Writer, this};
/// Handle for the output links between reco and gen particles
mutable DataHandle<edm4hep::TrackMCParticleLinkCollection> m_links{"TracksFromGenParticlesAlgAssociation", Gaudi::DataHandle::Writer, this};

// the next two functions have been copied from k4GaudiPandora and DDMarlinPandora / DDPandoraPFANewProcessor
// ideally they could be put in a separate file named e.g. DDUtils.cc
double getFieldFromCompact();
dd4hep::rec::LayeredCalorimeterData * getExtension(unsigned int includeFlag, unsigned int excludeFlag=0);
};

double TracksFromGenParticlesWithECalExtrapAlg::getFieldFromCompact() {
dd4hep::Detector& mainDetector = dd4hep::Detector::getInstance();
const double position[3]={0,0,0}; // position to calculate magnetic field at (the origin in this case)
double magneticFieldVector[3]={0,0,0}; // initialise object to hold magnetic field
mainDetector.field().magneticField(position,magneticFieldVector); // get the magnetic field vector from DD4hep
return magneticFieldVector[2]/dd4hep::tesla; // z component at (0,0,0)
}

dd4hep::rec::LayeredCalorimeterData * TracksFromGenParticlesWithECalExtrapAlg::getExtension(unsigned int includeFlag, unsigned int excludeFlag) {

dd4hep::rec::LayeredCalorimeterData * theExtension = 0;

dd4hep::Detector & mainDetector = dd4hep::Detector::getInstance();
const std::vector< dd4hep::DetElement>& theDetectors = dd4hep::DetectorSelector(mainDetector).detectors( includeFlag, excludeFlag );

debug() << " getExtension : includeFlag: " << dd4hep::DetType( includeFlag ) << " excludeFlag: " << dd4hep::DetType( excludeFlag )
<< " found : " << theDetectors.size() << " - first det: " << theDetectors.at(0).name()
<< endmsg;

if( theDetectors.size() != 1 ){

std::stringstream es ;
es << " getExtension: selection is not unique (or empty) includeFlag: " << dd4hep::DetType( includeFlag ) << " excludeFlag: " << dd4hep::DetType( excludeFlag )
<< " --- found detectors : " ;
for( unsigned i=0, N= theDetectors.size(); i<N ; ++i ){
es << theDetectors.at(i).name() << ", " ;
}
throw std::runtime_error( es.str() ) ;
}

theExtension = theDetectors.at(0).extension<dd4hep::rec::LayeredCalorimeterData>();

return theExtension;
}

TracksFromGenParticlesWithECalExtrapAlg::TracksFromGenParticlesWithECalExtrapAlg(const std::string& name, ISvcLocator* svcLoc) :
Gaudi::Algorithm(name, svcLoc) {
declareProperty("InputGenParticles", m_inputMCParticles, "input MCParticles");
declareProperty("OutputTracks", m_tracks, "Output tracks");
declareProperty("OutputMCRecoTrackParticleAssociation", m_links, "MCRecoTrackParticleAssociation");
m_Bz = 0.;
}

StatusCode TracksFromGenParticlesWithECalExtrapAlg::initialize() {
StatusCode sc = Gaudi::Algorithm::initialize();
if (sc.isFailure()) return sc;

// FIXME: handle exceptions if collections not found
for ( const auto& col : m_inputSimTrackerHitCollections ) {
debug() << "Creating handle for input SimTrackerHit collection : " << col << endmsg;
m_inputSimTrackerHitCollectionHandles.push_back(new DataHandle<edm4hep::SimTrackerHitCollection>(col, Gaudi::DataHandle::Reader, this));
}

// retrieve B field
m_Bz = getFieldFromCompact();
debug() << "B field (T) is : " << m_Bz << endmsg;

// retrieve tracker dimensions
//const dd4hep::rec::LayeredCalorimeterData * wrapperBarrelExtension = getExtension( ( dd4hep::DetType::TRACKER | dd4hep::DetType::BARREL ) , dd4hep::DetType::VERTEX );
//const dd4hep::rec::LayeredCalorimeterData * wrapperEndcapExtension = getExtension( ( dd4hep::DetType::TRACKER | dd4hep::DetType::ENDCAP ), dd4hep::DetType::VERTEX );
//const dd4hep::rec::LayeredCalorimeterData * vertexBarrelExtension = getExtension( ( dd4hep::DetType::TRACKER | dd4hep::DetType::BARREL | dd4hep::DetType::VERTEX ) );
//const dd4hep::rec::LayeredCalorimeterData * vertexEndcapExtension = getExtension( ( dd4hep::DetType::TRACKER | dd4hep::DetType::ENDCAP | dd4hep::DetType::VERTEX ) );

// retrieve ecal dimensions:
// - barrel: inner R, zmax
// - endcap: inner R, zmin, zmax
try {
const dd4hep::rec::LayeredCalorimeterData * eCalBarrelExtension = getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL),
( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) );
m_eCalBarrelInnerR = eCalBarrelExtension->extent[0] / dd4hep::mm;
m_eCalBarrelMaxZ = eCalBarrelExtension->extent[3] / dd4hep::mm;
debug() << "ECAL barrel extent: Rmin [mm] = " << m_eCalBarrelInnerR << endmsg;
debug() << "ECAL barrel extent: Zmax [mm] = " << m_eCalBarrelMaxZ << endmsg;
}
catch(...) {
warning() << "ECAL barrel extension not found" << endmsg;
m_eCalBarrelInnerR = 0.; // set to 0, will use it later to avoid projecting to the barrel
};

try {
const dd4hep::rec::LayeredCalorimeterData * eCalEndCapExtension = getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::ENDCAP),
( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) );
m_eCalEndCapInnerR = eCalEndCapExtension->extent[0] / dd4hep::mm;
m_eCalEndCapOuterR = eCalEndCapExtension->extent[1] / dd4hep::mm;
m_eCalEndCapInnerZ = eCalEndCapExtension->extent[2] / dd4hep::mm;
m_eCalEndCapOuterZ = eCalEndCapExtension->extent[3] / dd4hep::mm;
debug() << "ECAL endcap extent: Rmin [mm] = " << m_eCalEndCapInnerR << endmsg;
debug() << "ECAL endcap extent: Rmax [mm] = " << m_eCalEndCapOuterR << endmsg;
debug() << "ECAL endcap extent: Zmin [mm] = " << m_eCalEndCapInnerZ << endmsg;
debug() << "ECAL endcap extent: Zmax [mm] = " << m_eCalEndCapOuterZ << endmsg;
}
catch(...) {
warning() << "ECAL endcap extension not found" << endmsg;
m_eCalEndCapInnerR = 0.; // set to 0, will use it later to avoid projecting to the endcap
};
return StatusCode::SUCCESS;
}

StatusCode TracksFromGenParticlesWithECalExtrapAlg::execute(const EventContext&) const {
// Get the input MC particle collection
const edm4hep::MCParticleCollection* genParticleColl = m_inputMCParticles.get();

// Create the output track collection and track gen<->reco links collection
auto outputTrackCollection = new edm4hep::TrackCollection();
auto MCRecoTrackParticleAssociationCollection = new edm4hep::TrackMCParticleLinkCollection();

// loop over the gen particles, find charged ones, and create the corresponding reco particles
int iparticle = 0;
for (const auto& genParticle : *genParticleColl) {
debug() << endmsg;
debug() << "Gen. particle: " << genParticle << endmsg;
debug() <<" particle "<<iparticle++<<" PDG: "<< genParticle.getPDG() << " energy: "<<genParticle.getEnergy()
<< " charge: "<< genParticle.getCharge() << endmsg;
debug() << "Particle decayed in tracker: " << genParticle.isDecayedInTracker() << endmsg;

// consider only charged particles
if (genParticle.getCharge() == 0) continue;

// Building an helix out of MCParticle properties and B field
auto helixFromGenParticle = HelixClass_double();
auto vertex = genParticle.getVertex();
auto endpoint = genParticle.getEndpoint();
debug() << "Vertex radius: " << sqrt(vertex.x*vertex.x+vertex.y*vertex.y) << endmsg;
debug() << "Endpoint radius: " << sqrt(endpoint.x*endpoint.x+endpoint.y*endpoint.y) << endmsg;
double genParticleVertex[] = {vertex.x, vertex.y, vertex.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 at IP properties
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();
trackState_IP.referencePoint = edm4hep::Vector3f((float)genParticleVertex[0],(float)genParticleVertex[1],(float)genParticleVertex[2]);
trackFromGen.addToTrackStates(trackState_IP);

// find SimTrackerHits associated to genParticle, store hit position, momentum and time
std::vector<std::array<double,7> > trackHits;
for ( size_t ih=0; ih<m_inputSimTrackerHitCollectionHandles.size(); ih++ ) {
auto handle = dynamic_cast<DataHandle<edm4hep::SimTrackerHitCollection>*> (m_inputSimTrackerHitCollectionHandles[ih]);
const edm4hep::SimTrackerHitCollection* coll = handle->get();
for (const auto& hit : *coll) {
const edm4hep::MCParticle particle = hit.getParticle();
std::array<double,7> ahit{hit.x(), hit.y(), hit.z(), hit.getMomentum()[0], hit.getMomentum()[1], hit.getMomentum()[2], hit.getTime()};
if(particle.getObjectID() == genParticle.getObjectID()) trackHits.push_back(ahit);
}
}

if(!trackHits.empty())
{
// particles with at least one SimTrackerHit
debug() << "Number of SimTrackerHits: " << trackHits.size() << endmsg;

// sort the hits according to their time
std::sort(trackHits.begin(), trackHits.end(), [](const std::array<double,7> a, const std::array<double,7> b) {
return a[6]<b[6];
});

// TrackState at First Hit
auto trackState_AtFirstHit = edm4hep::TrackState {};
auto firstHit = trackHits.front();
double posAtFirstHit[] = {firstHit[0], firstHit[1], firstHit[2]};
double momAtFirstHit[] = {firstHit[3], firstHit[4], firstHit[5]};
debug() << "Radius of first hit: " << std::sqrt(firstHit[0]*firstHit[0] + firstHit[1]*firstHit[1]) << endmsg;
// get extrapolated momentum from the helix with ref point at IP
helixFromGenParticle.getExtrapolatedMomentum(posAtFirstHit,momAtFirstHit);
// produce new helix at first hit position
auto helixAtFirstHit = HelixClass_double();
helixAtFirstHit.Initialize_VP(posAtFirstHit, momAtFirstHit, genParticle.getCharge(), m_Bz);
// fill the TrackState parameters
trackState_AtFirstHit.location = edm4hep::TrackState::AtFirstHit;
trackState_AtFirstHit.D0 = helixAtFirstHit.getD0();
trackState_AtFirstHit.phi = helixAtFirstHit.getPhi0();
trackState_AtFirstHit.omega = helixAtFirstHit.getOmega();
trackState_AtFirstHit.Z0 = helixAtFirstHit.getZ0();
trackState_AtFirstHit.tanLambda = helixAtFirstHit.getTanLambda();
trackState_AtFirstHit.referencePoint = edm4hep::Vector3f((float)posAtFirstHit[0],(float)posAtFirstHit[1],(float)posAtFirstHit[2]);
trackFromGen.addToTrackStates(trackState_AtFirstHit);

// TrackState at Last Hit
auto trackState_AtLastHit = edm4hep::TrackState{};
auto lastHit = trackHits.back();
double posAtLastHit[] = {lastHit[0], lastHit[1], lastHit[2]};
double momAtLastHit[] = {lastHit[3], lastHit[4], lastHit[5]};
debug() << "Radius of last hit: " << std::sqrt(lastHit[0]*lastHit[0] + lastHit[1]*lastHit[1]) << endmsg;
// get extrapolated momentum from the helix with ref point at first hit
helixAtFirstHit.getExtrapolatedMomentum(posAtLastHit, momAtLastHit);
// produce new helix at last hit position
auto helixAtLastHit = HelixClass_double();
helixAtLastHit.Initialize_VP(posAtLastHit, momAtLastHit, genParticle.getCharge(), m_Bz);
// fill the TrackState parameters
trackState_AtLastHit.location = edm4hep::TrackState::AtLastHit;
trackState_AtLastHit.D0 = helixAtLastHit.getD0();
trackState_AtLastHit.phi = helixAtLastHit.getPhi0();
trackState_AtLastHit.omega = helixAtLastHit.getOmega();
trackState_AtLastHit.Z0 = helixAtLastHit.getZ0();
trackState_AtLastHit.tanLambda = helixAtLastHit.getTanLambda();
trackState_AtLastHit.referencePoint = edm4hep::Vector3f((float)posAtLastHit[0],
(float)posAtLastHit[1],
(float)posAtLastHit[2]);
// attach the TrackState to the track
trackFromGen.addToTrackStates(trackState_AtLastHit);

// TrackState at Calorimeter
if (m_eCalBarrelInnerR>0. || m_eCalEndCapInnerR>0.) {
auto trackState_AtCalorimeter = edm4hep::TrackState{};
pandora::CartesianVector bestECalProjection(0.f, 0.f, 0.f);

// create helix to project
const pandora::Helix helix(trackState_IP.phi,
trackState_IP.D0,
trackState_IP.Z0,
trackState_IP.omega,
trackState_IP.tanLambda,
m_Bz);
const pandora::CartesianVector& referencePoint(helix.GetReferencePoint());
const int signPz((helix.GetMomentum().GetZ() > 0.f) ? 1 : -1);

// First project to endcap
float minGenericTime(std::numeric_limits<float>::max());
if (m_eCalEndCapInnerR>0) {
(void)helix.GetPointInZ(static_cast<float>(signPz) * m_eCalEndCapInnerZ, referencePoint,
bestECalProjection, minGenericTime);
// GM: if the radius of the point on the plane corresponding to the endcap inner face is
// lower than the inner radius of the calorimeter, we might want to ignore it
// for the moment let's keep it as it might be useful for debugging the reconstruction
}

// Then project to barrel surface(s), and keep projection with lower arrival time
if (m_eCalBarrelInnerR>0) {
pandora::CartesianVector barrelProjection(0.f, 0.f, 0.f);
float genericTime(std::numeric_limits<float>::max());
const pandora::StatusCode statusCode(helix.GetPointOnCircle(m_eCalBarrelInnerR, referencePoint, barrelProjection, genericTime));
if ((pandora::STATUS_CODE_SUCCESS == statusCode) && (genericTime < minGenericTime)) {
minGenericTime = genericTime;
bestECalProjection = barrelProjection;
}
// GM: again, if the Z of the point on the cylinder of the barrel is beyond the
// max/min z of the detector, we might want to ignore it - but let's keep it
// for the moment as it might be useful for debugging the reconstruction
}

// get extrapolated position
double posAtCalorimeter[] = {bestECalProjection.GetX(), bestECalProjection.GetY(), bestECalProjection.GetZ()};
debug() << "Radius at calorimeter: " << std::sqrt(posAtCalorimeter[0]*posAtCalorimeter[0] + posAtCalorimeter[1]*posAtCalorimeter[1]) << endmsg;

// get extrapolated momentum from the helix with ref point at last hit
double momAtCalorimeter[] = {0.,0.,0.};
helixAtLastHit.getExtrapolatedMomentum(posAtCalorimeter, momAtCalorimeter);

// produce new helix at calorimeter position
auto helixAtCalorimeter = HelixClass_double();
helixAtCalorimeter.Initialize_VP(posAtCalorimeter, momAtCalorimeter, genParticle.getCharge(), m_Bz);

// fill the TrackState parameters
trackState_AtCalorimeter.location = edm4hep::TrackState::AtCalorimeter;
trackState_AtCalorimeter.D0 = helixAtCalorimeter.getD0();
trackState_AtCalorimeter.phi = helixAtCalorimeter.getPhi0();
trackState_AtCalorimeter.omega = helixAtCalorimeter.getOmega();
trackState_AtCalorimeter.Z0 = helixAtCalorimeter.getZ0();
trackState_AtCalorimeter.tanLambda = helixAtCalorimeter.getTanLambda();
trackState_AtCalorimeter.referencePoint = edm4hep::Vector3f((float)posAtCalorimeter[0],
(float)posAtCalorimeter[1],
(float)posAtCalorimeter[2]);
// attach the TrackState to the track
trackFromGen.addToTrackStates(trackState_AtCalorimeter);
}


outputTrackCollection->push_back(trackFromGen);

// Building the association between tracks and genParticles
auto MCRecoTrackParticleAssociation = edm4hep::MutableTrackMCParticleLink();
MCRecoTrackParticleAssociation.setFrom(trackFromGen);
MCRecoTrackParticleAssociation.setTo(genParticle);
MCRecoTrackParticleAssociationCollection->push_back(MCRecoTrackParticleAssociation);
}
}

// push the outputTrackCollection to event store
m_tracks.put(outputTrackCollection);
m_links.put(MCRecoTrackParticleAssociationCollection);

debug() << "Output tracks collection size: " << outputTrackCollection->size() << endmsg;

return StatusCode::SUCCESS;
}

StatusCode TracksFromGenParticlesWithECalExtrapAlg::finalize() { return Gaudi::Algorithm::finalize(); }


DECLARE_COMPONENT(TracksFromGenParticlesWithECalExtrapAlg)
4 changes: 2 additions & 2 deletions Tracking/test/runTracksFromGenParticles.py
Original file line number Diff line number Diff line change
@@ -3,7 +3,7 @@
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")
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/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml")

# Loading the output of the SIM step
from k4FWCore import IOSvc
@@ -24,7 +24,7 @@
from Configurables import PlotTrackHitDistances, RootHistSvc
from Configurables import Gaudi__Histograming__Sink__Root as RootHistoSink
plotTrackHitDistances = PlotTrackHitDistances("PlotTrackHitDistances",
InputSimTrackerHits = ["CDCHHits"],
InputSimTrackerHits = ["DCHCollection"],
InputTracksFromGenParticlesAssociation = tracksFromGenParticles.OutputMCRecoTrackParticleAssociation,
Bz = 2.0)

67 changes: 67 additions & 0 deletions Tracking/test/runTracksFromGenParticlesAlg.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
from k4FWCore import ApplicationMgr
from Gaudi.Configuration import INFO, 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 mu- --gun.thetaMin 20 --gun.thetaMax 160 --numberOfEvents 100 --outputFile ddsim_output_edm4hep.root --random.enableEventSeed --random.seed 42 --compactFile $K4GEO/FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.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"

# Geometry service
from Configurables import GeoSvc
import os
geoservice = GeoSvc("GeoSvc")
path_to_detector = os.environ.get("K4GEO", "")
detectors_to_use = [
'FCCee/ALLEGRO/compact/ALLEGRO_o1_v03/ALLEGRO_o1_v03.xml'
]
geoservice.detectors = [
os.path.join(path_to_detector, _det) for _det in detectors_to_use
]
geoservice.OutputLevel = INFO

# Calling TracksFromGenParticles
from Configurables import TracksFromGenParticlesWithECalExtrapAlg
tracksFromGenParticles = TracksFromGenParticlesWithECalExtrapAlg("TracksFromGenParticles",
InputGenParticles="MCParticles",
InputSimTrackerHits=[# "VertexBarrelCollection",
# "VertexEndcapCollection",
"DCHCollection",
"SiWrBCollection",
"SiWrDCollection"],
OutputTracks="TracksFromGenParticles",
OutputMCRecoTrackParticleAssociation="TracksFromGenParticlesAssociation",
OutputLevel=DEBUG)

# 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=["DCHCollection"],
InputTracksFromGenParticlesAssociation=["TracksFromGenParticlesAssociation"],
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, geoservice],
StopOnSignal=True,
)