From 5f27a439124414b924ec4d057264a99a9ba02150 Mon Sep 17 00:00:00 2001 From: jroloff Date: Tue, 9 Jul 2024 07:59:36 -0400 Subject: [PATCH] Adding tools for jet clustering for cluster and truth particles (#80) * Small changes to cmake lists * Adding tools for jet reconstruction * Adding truth jet algorithms * Fixing CMake change * Addressing some comments from the PR * Switching to thread safe Gaudi functional * Minor code cleanup * More code cleanup * Adding an extra algorithm to filter truth particles * One last removal for filtering * Adding separate jet clustering helper class * Adding some checks of the configuration in ClusterJet * A few updates to namespace etc * Addressing some changes from MR * Removing some commented lines * Update RecCalorimeter/src/components/FilterTruthParticles.cpp Co-authored-by: Brieuc Francois * Removing more unnecessary code * Adding tool descriptions * Addressing a few more MR comments * Fixing the tests in the CMake file * Changing the naming for truth particle filter * Fixing naming for truth jet test * Switching to using associations for truth jets --------- Co-authored-by: Brieuc Francois --- CMakeLists.txt | 4 +- RecCalorimeter/CMakeLists.txt | 22 +- .../include/RecCalorimeter/ClusterJet.h | 77 +++++ RecCalorimeter/src/components/ClusterJet.cpp | 47 +++ .../src/components/CreateCaloJet.cpp | 78 +++++ RecCalorimeter/src/components/CreateCaloJet.h | 59 ++++ .../src/components/CreateTruthJet.cpp | 67 ++++ .../src/components/CreateTruthJet.h | 59 ++++ .../FilterTruthParticlesForGenJets.cpp | 46 +++ .../FilterTruthParticlesForGenJets.h | 42 +++ RecCalorimeter/tests/options/runJetReco.py | 313 ++++++++++++++++++ .../tests/options/runTruthJetReco.py | 119 +++++++ cmake/FindFastJet.cmake | 29 ++ 13 files changed, 957 insertions(+), 5 deletions(-) create mode 100644 RecCalorimeter/include/RecCalorimeter/ClusterJet.h create mode 100644 RecCalorimeter/src/components/ClusterJet.cpp create mode 100644 RecCalorimeter/src/components/CreateCaloJet.cpp create mode 100644 RecCalorimeter/src/components/CreateCaloJet.h create mode 100644 RecCalorimeter/src/components/CreateTruthJet.cpp create mode 100644 RecCalorimeter/src/components/CreateTruthJet.h create mode 100644 RecCalorimeter/src/components/FilterTruthParticlesForGenJets.cpp create mode 100644 RecCalorimeter/src/components/FilterTruthParticlesForGenJets.h create mode 100644 RecCalorimeter/tests/options/runJetReco.py create mode 100644 RecCalorimeter/tests/options/runTruthJetReco.py create mode 100644 cmake/FindFastJet.cmake diff --git a/CMakeLists.txt b/CMakeLists.txt index e2bd08ca..86033496 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,6 +29,7 @@ find_package(EDM4HEP 0.10.1) # implicit: Podio find_package(DD4hep) find_package(k4geo) find_package(ONNXRuntime) +find_package(FastJet) # New versions of ONNRuntime package provide onnxruntime-Config.cmake # and use the name onnxruntime find_package(onnxruntime) @@ -64,9 +65,6 @@ add_subdirectory(RecFCChhCalorimeter) add_subdirectory(RecFCCeeCalorimeter) - - - install(EXPORT ${PROJECT_NAME}Targets NAMESPACE ${PROJECT_NAME}:: FILE "${PROJECT_NAME}Targets.cmake" diff --git a/RecCalorimeter/CMakeLists.txt b/RecCalorimeter/CMakeLists.txt index e1e117a7..70811457 100644 --- a/RecCalorimeter/CMakeLists.txt +++ b/RecCalorimeter/CMakeLists.txt @@ -3,9 +3,10 @@ ################################################################################ -file(GLOB _module_sources src/components/*.cpp) +file(GLOB _module_sources src/components/*.cpp ) + gaudi_add_module(k4RecCalorimeterPlugins - SOURCES ${_module_sources} + SOURCES ${_module_sources} LINK k4FWCore::k4FWCore k4FWCore::k4Interface Gaudi::GaudiAlgLib @@ -17,6 +18,7 @@ gaudi_add_module(k4RecCalorimeterPlugins ROOT::Hist k4geo::detectorSegmentations k4geo::detectorCommon + ${FASTJET_LIBRARIES} ) install(TARGETS k4RecCalorimeterPlugins @@ -26,8 +28,24 @@ install(TARGETS k4RecCalorimeterPlugins COMPONENT dev) + +include_directories(${FASTJET_INCLUDE_DIRS} include/RecCalorimeter) +link_directories(${FASTJET_LIBRARIES}) + include(CTest) +add_test(NAME FCC_createJets + COMMAND k4run RecCalorimeter/tests/options/runJetReco.py + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} +) +set_test_env(FCC_createJets) + +add_test(NAME FCC_createTruthJets + COMMAND k4run RecCalorimeter/tests/options/runTruthJetReco.py + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} +) +set_test_env(FCC_createTruthJets) + #install(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/tests/options DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}/Reconstruction/RecCalorimeter) # #gaudi_add_test(genJetClustering diff --git a/RecCalorimeter/include/RecCalorimeter/ClusterJet.h b/RecCalorimeter/include/RecCalorimeter/ClusterJet.h new file mode 100644 index 00000000..f2b35172 --- /dev/null +++ b/RecCalorimeter/include/RecCalorimeter/ClusterJet.h @@ -0,0 +1,77 @@ +#ifndef RECCALORIMETER_CLUSTERJET_H +#define RECCALORIMETER_CLUSTERJET_H + +// std +#include + +// Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" + +#include "fastjet/JetDefinition.hh" + +// EDM4hep +#include "edm4hep/ClusterCollection.h" + + +namespace k4::recCalo{ + +class ClusterInfo : public fastjet::PseudoJet::UserInfoBase{ + public: + ClusterInfo(const int & index) : _index(index){} + int index() const{ return _index; } + protected: + int _index; +}; + + +/** @class ClusterJet k4RecCalorimeter/RecCalorimeter/src/components/ClusterJet.h + * + * Helper class for reconstructing jets, used to avoid code duplication in the classes that require specific input types. + * It takes a vector of fastjet::PseudoJets as inputs, and returns the a vector of pseudojets after clustering. + * + * JetAlg: A string corresponding to the jet algorithm for clustering + * JetRadius: The radius of the jets being clustered + * MinPt: The pT threshold below which jets are ignored + * isExclusiveClustering: 1 if jets should use an exclusive clustering, 0 otherwise + * clusterArgs: Any other clustering arguments, such as the number of jets for exclusive clustering. + * This is not used everywhere, so check the code to make sure it is implemented for your use-case. + * + * @author Jennifer Roloff + * @date 2024-7 + */ + + + +class ClusterJet { +public: + ClusterJet(std::string jetAlg, double jetRadius, int isExclusiveClustering = 0, double minPt = 0, int clusterArgs = 0); + StatusCode initialize(); + ~ClusterJet(){if(m_clustSeq) delete m_clustSeq;} + + std::vector cluster(const std::vector clustersPJ); + +private: + + + + std::map m_jetAlgMap = { + {"kt", fastjet::JetAlgorithm::kt_algorithm}, + {"cambridge", fastjet::JetAlgorithm::cambridge_algorithm}, + {"antikt", fastjet::JetAlgorithm::antikt_algorithm}, + {"genkt", fastjet::JetAlgorithm::genkt_algorithm}, + {"ee_kt", fastjet::JetAlgorithm::ee_kt_algorithm}, + {"ee_genkt", fastjet::JetAlgorithm::ee_genkt_algorithm}, + }; + + std::string m_jetAlg = "antikt"; + double m_jetRadius = 0.4; + int m_isExclusiveClustering = 0; // Inclusive clustering by default + double m_minPt = 10; // Only relevant for inclusive clustering + int m_njets = 0; // Only relevant for exclusive clustering + + fastjet::ClusterSequence* m_clustSeq; + +}; + +} +#endif /* RECCALORIMETER_CLUSTERJET_H */ diff --git a/RecCalorimeter/src/components/ClusterJet.cpp b/RecCalorimeter/src/components/ClusterJet.cpp new file mode 100644 index 00000000..e1838337 --- /dev/null +++ b/RecCalorimeter/src/components/ClusterJet.cpp @@ -0,0 +1,47 @@ +#include "ClusterJet.h" + +// std +#include +#include + +namespace k4::recCalo{ + +ClusterJet::ClusterJet(std::string jetAlg, double jetRadius, int isExclusiveClustering, double minPt, int njets): m_jetAlg(jetAlg), m_jetRadius(jetRadius), m_isExclusiveClustering(isExclusiveClustering), m_minPt(minPt), m_njets(njets){ + m_clustSeq = nullptr; +} + + +StatusCode ClusterJet::initialize(){ + if (m_jetAlgMap.find(m_jetAlg) == m_jetAlgMap.end()) { + std::cout << "ERROR: " << " is not in the list of supported jet algorithms" << std::endl;; + return StatusCode::FAILURE; + } + if(m_isExclusiveClustering > 1){ + std::cout << "ERROR: " << "Clustering of " << m_isExclusiveClustering << " is currently not supported" << std::endl;; + return StatusCode::FAILURE; + } + return StatusCode::SUCCESS; +} + +std::vector ClusterJet::cluster(const std::vector clustersPJ) { + std::vector jets; + + fastjet::JetDefinition* jetDef = new fastjet::JetDefinition( m_jetAlgMap.at(m_jetAlg), m_jetRadius); + if(m_clustSeq) delete m_clustSeq; + m_clustSeq = new fastjet::ClusterSequence(clustersPJ, *jetDef); + + // Note: initialize has already checked if m_isExclusiveClustering has the right range + if(m_isExclusiveClustering == 0){ + jets = fastjet::sorted_by_pt(m_clustSeq->inclusive_jets(m_minPt)); + } + else if(m_isExclusiveClustering == 1){ + jets = fastjet::sorted_by_pt(m_clustSeq->exclusive_jets(m_njets)); + } + + return jets; + +} + + +} + diff --git a/RecCalorimeter/src/components/CreateCaloJet.cpp b/RecCalorimeter/src/components/CreateCaloJet.cpp new file mode 100644 index 00000000..3de6c81f --- /dev/null +++ b/RecCalorimeter/src/components/CreateCaloJet.cpp @@ -0,0 +1,78 @@ +#include "CreateCaloJet.h" + +// std +#include +#include + +// EDM4hep +#include "edm4hep/ClusterCollection.h" + +DECLARE_COMPONENT(CreateCaloJet) + +CreateCaloJet::CreateCaloJet(const std::string& name, ISvcLocator* svcLoc) : Transformer(name, svcLoc, + KeyValue("InputCollection", "CorrectedCaloClusters"), + KeyValue("OutputCollection", "Jets")) { + declareProperty("JetAlg", m_jetAlg, "Name of jet clustering algorithm"); + declareProperty("JetRadius", m_jetRadius, "Jet clustering radius"); + declareProperty("MinPt", m_minPt, "Minimum pT for saved jets"); + declareProperty("isExclusiveClustering", m_isExclusive, "1 if exclusive, 0 if inclusive"); + +} + + +StatusCode CreateCaloJet::initialize() { + clusterer = new k4::recCalo::ClusterJet(m_jetAlg, m_jetRadius, m_isExclusive, m_minPt); + + + return clusterer->initialize(); +} + +edm4hep::ReconstructedParticleCollection CreateCaloJet::operator()(const edm4hep::ClusterCollection& input) const{ + std::vector clustersPJ; + int i=0; + + // Need to convert cluster position in space to the momentum + for(auto cluster: input){ + const edm4hep::Vector3f position = cluster.getPosition(); + float x = position.x; + float y = position.y; + float z = position.z; + double theta = acos(sqrt(z*z / (x*x + y*y + z*z))); + double eta = -log(tan(theta / 2.)); + double phi = atan2(y, x); + double pT = cluster.getEnergy()* sqrt( (x*x + y*y) / (x*x + y*y + z*z)); + + // Take clusters to be massless + fastjet::PseudoJet clusterPJ; + clusterPJ.reset_momentum_PtYPhiM(pT, eta, phi, 0); + // Attach the cluster index to the pseudojet + clusterPJ.set_user_info(new k4::recCalo::ClusterInfo(i)); + clustersPJ.push_back(clusterPJ); + i++; + } + + + std::vector inclusiveJets = clusterer->cluster(clustersPJ); + + + edm4hep::ReconstructedParticleCollection edmJets = edm4hep::ReconstructedParticleCollection(); + //Add a reconstructed particle for each jet + for(auto cjet : inclusiveJets){ + edm4hep::MutableReconstructedParticle jet; + jet.setMomentum(edm4hep::Vector3f(cjet.px(), cjet.py(), cjet.pz())); + jet.setEnergy(cjet.e()); + jet.setMass(cjet.m()); + + // Also add the clusters that were used to make the jets + std::vector constits = cjet.constituents(); + for(auto constit : constits){ + int index = constit.user_info().index(); + jet.addToClusters((input)[index]); + } + edmJets.push_back(jet); + } + + return edmJets; +} + + diff --git a/RecCalorimeter/src/components/CreateCaloJet.h b/RecCalorimeter/src/components/CreateCaloJet.h new file mode 100644 index 00000000..94f4265a --- /dev/null +++ b/RecCalorimeter/src/components/CreateCaloJet.h @@ -0,0 +1,59 @@ +#ifndef RECCALORIMETER_CREATECALOJET_H +#define RECCALORIMETER_CREATECALOJET_H + + +// Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiAlg/Transformer.h" +#include "k4FWCore/BaseClass.h" + +// k4FWCore +#include "k4FWCore/DataHandle.h" + +//EDM4Hep +#include "edm4hep/ReconstructedParticleCollection.h" + +// EDM4hep +#include "edm4hep/ClusterCollection.h" + + +#include "ClusterJet.h" + + + +/** @class CreateCaloJet k4RecCalorimeter/RecFCCeeCalorimeter/src/components/CreateCaloJet.h + * + * Tool for reconstructing jets from calorimeter clusters. + * It takes as input a ClusterCollection, and it outputs a ReconstructedParticleCollection with the jets. + * + * JetAlg: A string corresponding to the jet algorithm for clustering + * JetRadius: The radius of the jets being clustered + * MinPt: The pT threshold below which jets are ignored + * isExclusiveClustering: 1 if jets should use an exclusive clustering, 0 otherwise + * + * @author Jennifer Roloff + * @date 2024-7 + */ + + + +class CreateCaloJet : public Gaudi::Functional::Transformer { +public: + CreateCaloJet(const std::string& name, ISvcLocator* svcLoc); + edm4hep::ReconstructedParticleCollection operator()(const edm4hep::ClusterCollection& input) const override; + + StatusCode initialize() override; + +private: + + + double m_minPt = 10; + double m_jetRadius = 0.4; + std::string m_jetAlg = "antikt"; + int m_isExclusive = 0; + + k4::recCalo::ClusterJet* clusterer; + + +}; +#endif /* RECCALORIMETER_CREATECALOJET_H */ diff --git a/RecCalorimeter/src/components/CreateTruthJet.cpp b/RecCalorimeter/src/components/CreateTruthJet.cpp new file mode 100644 index 00000000..f2ea535e --- /dev/null +++ b/RecCalorimeter/src/components/CreateTruthJet.cpp @@ -0,0 +1,67 @@ +#include "CreateTruthJet.h" + +// std +#include +#include + +DECLARE_COMPONENT(CreateTruthJet) + +CreateTruthJet::CreateTruthJet(const std::string& name, ISvcLocator* svcLoc) : Transformer(name, svcLoc, + KeyValue("InputCollection", "MCParticles"), + KeyValue("OutputCollection", "TruthJets")) { + declareProperty("JetAlg", m_jetAlg, "Name of jet clustering algorithm"); + declareProperty("JetRadius", m_jetRadius, "Jet clustering radius"); + declareProperty("MinPt", m_minPt, "Minimum pT for saved jets"); + declareProperty("isExclusiveClustering", m_isExclusive, "1 if exclusive, 0 if inclusive"); + declareProperty("outputAssociation", m_outputAssocCollection, "TruthJetParticleAssociation"); +} + + + +StatusCode CreateTruthJet::initialize() { + clusterer = new k4::recCalo::ClusterJet(m_jetAlg, m_jetRadius, m_isExclusive, m_minPt); + return clusterer->initialize(); + +} + +edm4hep::ReconstructedParticleCollection CreateTruthJet::operator()(const edm4hep::MCParticleCollection& input) const{ + std::vector clustersPJ; + auto& assoc = *(m_outputAssocCollection.createAndPut()); + int i=0; + for(auto particle: input){ + fastjet::PseudoJet clusterPJ(particle.getMomentum().x, particle.getMomentum().y, particle.getMomentum().z, particle.getEnergy()); + clusterPJ.set_user_info(new k4::recCalo::ClusterInfo(i)); + clustersPJ.push_back(clusterPJ); + i++; + } + + std::vector inclusiveJets = clusterer->cluster(clustersPJ); + + + edm4hep::ReconstructedParticleCollection edmJets = edm4hep::ReconstructedParticleCollection(); + + for(auto cjet : inclusiveJets){ + edm4hep::MutableReconstructedParticle jet; + jet.setMomentum(edm4hep::Vector3f(cjet.px(), cjet.py(), cjet.pz())); + jet.setEnergy(cjet.e()); + jet.setMass(cjet.m()); + + std::vector constits = cjet.constituents(); + + for(auto constit : constits){ + int index = constit.user_info().index(); + + edm4hep::MutableMCRecoParticleAssociation association; + association.setRec(jet); + association.setSim((input)[index]); + assoc.push_back(association); + } + + edmJets.push_back(jet); + } + + return edmJets; +} + + + diff --git a/RecCalorimeter/src/components/CreateTruthJet.h b/RecCalorimeter/src/components/CreateTruthJet.h new file mode 100644 index 00000000..cce2362e --- /dev/null +++ b/RecCalorimeter/src/components/CreateTruthJet.h @@ -0,0 +1,59 @@ +#ifndef RECCALORIMETER_CREATETRUTHJET_H +#define RECCALORIMETER_CREATETRUTHJET_H + +// Gaudi +#include "GaudiKernel/ToolHandle.h" +#include "GaudiAlg/Transformer.h" +#include "k4FWCore/BaseClass.h" + +// k4FWCore +#include "k4FWCore/DataHandle.h" + +#include "edm4hep/ReconstructedParticleCollection.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/MCRecoParticleAssociation.h" +#include "edm4hep/MCRecoParticleAssociationCollection.h" + +#include "ClusterJet.h" + + +/** @class CreateTruthJet k4RecCalorimeter/RecCalorimeter/src/components/CreateTruthJet.h + * + * Tool for reconstructing jets from truth particles. + * It takes as input an MCParticleCollection, and it outputs a ReconstructedParticleCollection with the jets. + * + * JetAlg: A string corresponding to the jet algorithm for clustering + * JetRadius: The radius of the jets being clustered + * MinPt: The pT threshold below which jets are ignored + * isExclusiveClustering: 1 if jets should use an exclusive clustering, 0 otherwise + * outputAssociation: Name of the output association collection to link jets to their constituents + * + + * @author Jennifer Roloff + * @date 2024-7 + */ + + +class CreateTruthJet : public Gaudi::Functional::Transformer { +public: + CreateTruthJet(const std::string& name, ISvcLocator* svcLoc); + edm4hep::ReconstructedParticleCollection operator()(const edm4hep::MCParticleCollection& input) const override final; + + StatusCode initialize() override final; + +private: + + + mutable DataHandle m_outputAssocCollection{"MCRecoParticleAssociation", + Gaudi::DataHandle::Writer, this}; + + double m_minPt = 1; + double m_jetRadius = 0.4; + std::string m_jetAlg = "antikt"; + int m_isExclusive = 0; + + k4::recCalo::ClusterJet* clusterer; + + +}; +#endif /* RECCALORIMETER_CREATETRUTHJET_H */ diff --git a/RecCalorimeter/src/components/FilterTruthParticlesForGenJets.cpp b/RecCalorimeter/src/components/FilterTruthParticlesForGenJets.cpp new file mode 100644 index 00000000..5bb899d5 --- /dev/null +++ b/RecCalorimeter/src/components/FilterTruthParticlesForGenJets.cpp @@ -0,0 +1,46 @@ +#include "FilterTruthParticlesForGenJets.h" + +// std +#include +#include + + +DECLARE_COMPONENT(FilterTruthParticlesForGenJets) + +FilterTruthParticlesForGenJets::FilterTruthParticlesForGenJets(const std::string& name, ISvcLocator* svcLoc) : Transformer(name, svcLoc, + KeyValue("InputCollection", "MCParticles"), + KeyValue("OutputCollection", "FilteredMCParticles")) { + + declareProperty("SaveMuons", m_useMuons, "Bool for if muons are stored"); + declareProperty("MinPt", m_minPt, "Minimum pT for filtered particles"); +} + + + +edm4hep::MCParticleCollection FilterTruthParticlesForGenJets::operator()(const edm4hep::MCParticleCollection& input) const{ + + edm4hep::MCParticleCollection filteredParticles = edm4hep::MCParticleCollection(); + + for(auto particle: input){ + + // Only want to include final state particles (before GEANT simulation) + if(particle.isCreatedInSimulation()) continue; + if(particle.getGeneratorStatus() != 1) continue; + + // Sometimes we want to ignore muons for jet reconstruction + if(!m_useMuons && std::abs(particle.getPDG()) == 13) continue; + + // We may want to impose a minimum pT requirement + double input_particle_pt = std::sqrt(particle.getMomentum().x * particle.getMomentum().x + particle.getMomentum().y * particle.getMomentum().y); + if(input_particle_pt < m_minPt) continue; + + edm4hep::MutableMCParticle filteredParticle = particle.clone(); + filteredParticles.push_back(filteredParticle); + } + + + return filteredParticles; +} + + + diff --git a/RecCalorimeter/src/components/FilterTruthParticlesForGenJets.h b/RecCalorimeter/src/components/FilterTruthParticlesForGenJets.h new file mode 100644 index 00000000..ef75ffbb --- /dev/null +++ b/RecCalorimeter/src/components/FilterTruthParticlesForGenJets.h @@ -0,0 +1,42 @@ +#ifndef RECCALORIMETER_FILTERTRUTHPARTICLESFORGENJETS_H +#define RECCALORIMETER_FILTERTRUTHPARTICLESFORGENJETS_H + + +// Gaudi +#include "GaudiAlg/Transformer.h" +#include "k4FWCore/BaseClass.h" + + +#include "edm4hep/MCParticleCollection.h" + +#include "fastjet/JetDefinition.hh" + + +/** @class FilterTruthParticlesForGenJets k4RecCalorimeter/RecCalorimeter/src/components/FilterTruthParticlesForGenJets.h + * + * Tool to filter truth particles based on their status code. + * Only detector stable particles are stored in a new MCParticleCollection. + * This can provide inputs to truth jet clustering. + * + * MinPt: the minimum pT of particles that are stored + * SaveMuons: 1 if muons should be included in the collection, 0 otherwise. + * + * @author Jennifer Roloff + * @date 2024-7 + */ + + +class FilterTruthParticlesForGenJets : public Gaudi::Functional::Transformer { +public: + FilterTruthParticlesForGenJets(const std::string& name, ISvcLocator* svcLoc); + edm4hep::MCParticleCollection operator()(const edm4hep::MCParticleCollection& input) const override; + + +private: + + double m_minPt = 0.0; + int m_useMuons = 0; + + +}; +#endif /* RECCALORIMETER_FILTERTRUTHPARTICLES_H */ diff --git a/RecCalorimeter/tests/options/runJetReco.py b/RecCalorimeter/tests/options/runJetReco.py new file mode 100644 index 00000000..3e6d6e57 --- /dev/null +++ b/RecCalorimeter/tests/options/runJetReco.py @@ -0,0 +1,313 @@ +import os +from Gaudi.Configuration import * +from GaudiKernel import SystemOfUnits as units +from GaudiKernel import PhysicalConstants as constants +from GaudiKernel.SystemOfUnits import MeV, GeV, tesla + +from Configurables import ApplicationMgr +ApplicationMgr().EvtSel = 'NONE' +ApplicationMgr().EvtMax = 2 +ApplicationMgr().OutputLevel = INFO +ApplicationMgr().StopOnSignal = True +ApplicationMgr().ExtSvc += ['RndmGenSvc'] + +from Configurables import k4DataSvc +## Data service +podioevent = k4DataSvc("EventDataSvc") +ApplicationMgr().ExtSvc += [podioevent] + +from Configurables import MomentumRangeParticleGun +guntool = MomentumRangeParticleGun() +guntool.ThetaMin = 45 * constants.pi / 180. +guntool.ThetaMax = 135 * constants.pi / 180. +guntool.PhiMin = 0. +guntool.PhiMax = 2. * constants.pi +guntool.MomentumMin = 1. *units.GeV +guntool.MomentumMax = 1. *units.GeV +guntool.PdgCodes = [11] # 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ + +from Configurables import GenAlg +gen = GenAlg() +gen.SignalProvider=guntool +gen.hepmc.Path = "hepmc" +ApplicationMgr().TopAlg += [gen] + +from Configurables import HepMCToEDMConverter +## reads an HepMC::GenEvent from the data service and writes a collection of EDM Particles +hepmc_converter = HepMCToEDMConverter("Converter") +hepmc_converter.hepmc.Path="hepmc" +hepmc_converter.GenParticles.Path="GenParticles" +ApplicationMgr().TopAlg += [hepmc_converter] + + + +################## Simulation setup +# Detector geometry +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +# if FCC_DETECTORS is empty, this should use relative path to working directory +path_to_detectors = os.environ.get("K4GEO", "") +detectors = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v01/ALLEGRO_o1_v01.xml' +] +# prefix all xmls with path_to_detectors +for det in detectors: + geoservice.detectors += [os.path.join(path_to_detectors, det)] +geoservice.OutputLevel = INFO +ApplicationMgr().ExtSvc += [geoservice] + + +# Geant4 service +# Configures the Geant simulation: geometry, physics list and user actions +from Configurables import SimG4Svc +# giving the names of tools will initialize the tools of that type +geantservice = SimG4Svc("SimG4Svc") +geantservice.detector = "SimG4DD4hepDetector" +geantservice.physicslist = "SimG4FtfpBert" +geantservice.actions = "SimG4FullSimActions" +# Fixed seed to have reproducible results, change it for each job if you split one production into several jobs +# Mind that if you leave Gaudi handle random seed and some job start within the same second (very likely) you will have duplicates +geantservice.randomNumbersFromGaudi = False +geantservice.seedValue = 4242 +# Range cut +geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] +ApplicationMgr().ExtSvc += [geantservice] + + + +# Magnetic field +from Configurables import SimG4ConstantMagneticFieldTool +field = SimG4ConstantMagneticFieldTool("SimG4ConstantMagneticFieldTool") +field.FieldComponentZ = -2 * units.tesla +field.FieldOn = True +field.IntegratorStepper="ClassicalRK4" + +# Geant4 algorithm +# Translates EDM to G4Event, passes the event to G4, writes out outputs via tools +# and a tool that saves the calorimeter hits + +# Detector readouts +# ECAL +ecalBarrelReadoutName = "ECalBarrelEta" +ecalBarrelReadoutNamePhiEta = "ECalBarrelPhiEta" +ecalEndcapReadoutName = "ECalEndcapPhiEta" +# HCAL +#hcalReadoutName = "HCalBarrelReadout" +#extHcalReadoutName = "HCalExtBarrelReadout" + +from Configurables import SimG4Alg +from Configurables import SimG4SaveCalHits +SimG4Alg("SimG4Alg").outputs = [] + +saveECalBarrelTool = SimG4SaveCalHits("saveECalBarrelHits") +saveECalBarrelTool.readoutNames = ["ECalBarrelEta"] +saveECalBarrelTool.CaloHits.Path = "ECalBarrelPositionedHits" +SimG4Alg("SimG4Alg").outputs += [saveECalBarrelTool] + +saveECalEndcapTool = SimG4SaveCalHits("saveECalEndcapHits") +saveECalEndcapTool.readoutNames = ["ECalEndcapPhiEta"] +saveECalEndcapTool.CaloHits.Path = "ECalEndcapPositionedHits" +SimG4Alg("SimG4Alg").outputs += [saveECalEndcapTool] + +#saveHCalTool = SimG4SaveCalHits("saveHCalBarrelHits") +#saveHCalTool.readoutNames = ["HCalBarrelReadout"] +#saveHCalTool.CaloHits.Path = "HCalBarrelPositionedHits" +#SimG4Alg("SimG4Alg").outputs += [saveHCalTool] + +# next, create the G4 algorithm, giving the list of names of tools ("XX/YY") +from Configurables import SimG4PrimariesFromEdmTool +particle_converter = SimG4PrimariesFromEdmTool("EdmConverter") +particle_converter.GenParticles.Path = "GenParticles" + +from Configurables import SimG4Alg +geantsim = SimG4Alg("SimG4Alg") +geantsim.eventProvider = particle_converter +ApplicationMgr().TopAlg += [geantsim] + +# Uncomment the following if history from Geant4 decays is needed (e.g. to get the photons from pi0) +#from Configurables import SimG4FullSimActions, SimG4SaveParticleHistory +#actions = SimG4FullSimActions() +#actions.enableHistory = True +#actions.energyCut = 0.2 * GeV +#saveHistTool = SimG4SaveParticleHistory("saveHistory") +#SimG4Alg("SimG4Alg").outputs += [saveHistTool] +#geantservice.actions = actions + + +############## Digitization (Merging hits into cells, EM scale calibration) +# EM scale calibration (sampling fraction) +from Configurables import CalibrateInLayersTool + +calibEcalBarrel = CalibrateInLayersTool("CalibrateECalBarrel") +calibEcalBarrel.samplingFraction = [0.3632447480841956] * 1 + [0.13187261040190248] * 1 + [0.14349714292943705] * 1 + [0.150266118277841] * 1 + [0.15502683375826457] * 1 + [0.15954408786354762] * 1 + [0.16375302347299436] * 1 + [0.16840384714588075] * 1 + [0.17219540619311383] * 1 + [0.1755068643940401] * 1 + [0.17816980262822366] * 1 + [0.18131266048670405] * 1 +calibEcalBarrel.readoutName = "ECalBarrelEta" +calibEcalBarrel.layerFieldName = "layer" + +from Configurables import CalibrateCaloHitsTool +#calibHcells = CalibrateCaloHitsTool("CalibrateHCal") +#calibHcells.invSamplingFraction = "41.66" +calibEcalEndcap = CalibrateCaloHitsTool("CalibrateECalEndcap") +calibEcalEndcap.invSamplingFraction="4.27" + +# Create cells in ECal barrel +# 1. step - merge hits into cells with Eta and module segmentation (phi module is a 'physical' cell i.e. lead + LAr + PCB + LAr +lead) +# 2. step - rewrite the cellId using the Eta-Phi segmentation (merging several modules into one phi readout cell) +from Configurables import CreateCaloCells +createEcalBarrelCellsStep1 = CreateCaloCells("CreateECalBarrelCellsStep1") +createEcalBarrelCellsStep1.doCellCalibration = True +createEcalBarrelCellsStep1.calibTool = calibEcalBarrel +createEcalBarrelCellsStep1.addCellNoise = False +createEcalBarrelCellsStep1.filterCellNoise = False +# todo: add when update on cvmfs +createEcalBarrelCellsStep1.addPosition = True +createEcalBarrelCellsStep1.hits = "ECalBarrelPositionedHits" +createEcalBarrelCellsStep1.cells = "ECalBarrelCellsStep1" +ApplicationMgr().TopAlg += [createEcalBarrelCellsStep1] + +## Use Phi-Theta segmentation in ECal barrel +from Configurables import RedoSegmentation +resegmentEcalBarrel = RedoSegmentation("ReSegmentationEcal") +# old bitfield (readout) +resegmentEcalBarrel.oldReadoutName = "ECalBarrelEta" +# specify which fields are going to be altered (deleted/rewritten) +resegmentEcalBarrel.oldSegmentationIds = ["module"] +# new bitfield (readout), with new segmentation +resegmentEcalBarrel.newReadoutName = "ECalBarrelPhiEta" +resegmentEcalBarrel.inhits = "ECalBarrelCellsStep1" +resegmentEcalBarrel.outhits = "ECalBarrelCellsStep2" +ApplicationMgr().TopAlg += [resegmentEcalBarrel] + +EcalBarrelCellsName = "ECalBarrelCells" +from Configurables import CreateCaloCells +createEcalBarrelCells = CreateCaloCells("CreateECalBarrelCells") +createEcalBarrelCells.doCellCalibration = False +createEcalBarrelCells.addCellNoise = False +createEcalBarrelCells.filterCellNoise = False +createEcalBarrelCells.hits = "ECalBarrelCellsStep2" +createEcalBarrelCells.cells = "ECalBarrelCells" +ApplicationMgr().TopAlg += [createEcalBarrelCells] + +# Ecal barrel cell positions (good for physics, all coordinates set properly) +from Configurables import CellPositionsECalBarrelTool +cellPositionEcalBarrelTool = CellPositionsECalBarrelTool("CellPositionsECalBarrel") +cellPositionEcalBarrelTool.readoutName = "ECalBarrelPhiEta" + +from Configurables import CreateCaloCellPositionsFCCee +createEcalBarrelPositionedCells = CreateCaloCellPositionsFCCee("ECalBarrelPositionedCells") +createEcalBarrelPositionedCells.positionsTool = cellPositionEcalBarrelTool +createEcalBarrelPositionedCells.hits.Path = "ECalBarrelCells" +createEcalBarrelPositionedCells.positionedHits.Path = "ECalBarrelPositionedCells" +ApplicationMgr().TopAlg += [createEcalBarrelPositionedCells] + +# Create cells in HCal +# 1. step - merge hits into cells with the default readout +#createHcalBarrelCells = CreateCaloCells("CreateHCaloCells") +#createHcalBarrelCells.doCellCalibration = True +#createHcalBarrelCells.calibTool = calibHcells +#createHcalBarrelCells.addPosition = True +#createHcalBarrelCells.addCellNoise = False +#createHcalBarrelCells.filterCellNoise = False +#createHcalBarrelCells.hits = "HCalBarrelPositionedHits" +#createHcalBarrelCells.cells = "HCalBarrelCells" +#ApplicationMgr().TopAlg += [createHcalBarrelCells] + +# Create cells in ECal Endcaps +createEcalEndcapCells = CreateCaloCells("CreateEcalEndcapCaloCells") +createEcalEndcapCells.doCellCalibration=True +createEcalEndcapCells.calibTool=calibEcalEndcap +createEcalEndcapCells.addCellNoise=False +createEcalEndcapCells.filterCellNoise=False +createEcalEndcapCells.hits.Path="ECalEndcapPositionedHits" +createEcalEndcapCells.cells.Path="ECalEndcapCells" +ApplicationMgr().TopAlg += [createEcalEndcapCells] + +#Empty cells for parts of calorimeter not implemented yet +from Configurables import CreateEmptyCaloCellsCollection +createemptycells = CreateEmptyCaloCellsCollection("CreateEmptyCaloCells") +createemptycells.cells.Path = "emptyCaloCells" +ApplicationMgr().TopAlg += [createemptycells] + +from Configurables import CaloTowerTool +towers = CaloTowerTool("towers") +towers.deltaEtaTower = 0.01 +towers.deltaPhiTower = 2 * constants.pi / 768. +towers.ecalBarrelReadoutName = ecalBarrelReadoutNamePhiEta +towers.ecalEndcapReadoutName = "" +towers.ecalFwdReadoutName = "" +towers.hcalBarrelReadoutName = "" +towers.hcalExtBarrelReadoutName = "" +towers.hcalEndcapReadoutName = "" +towers.hcalFwdReadoutName = "" +towers.ecalBarrelCells.Path = EcalBarrelCellsName +towers.ecalEndcapCells.Path = "emptyCaloCells" +towers.ecalFwdCells.Path = "emptyCaloCells" +towers.hcalBarrelCells.Path = "emptyCaloCells" +towers.hcalExtBarrelCells.Path = "emptyCaloCells" +towers.hcalEndcapCells.Path = "emptyCaloCells" +towers.hcalFwdCells.Path = "emptyCaloCells" + + + +from Configurables import CreateCaloClustersSlidingWindow +createClusters = CreateCaloClustersSlidingWindow("CreateClusters") +# Cluster variables +createClusters.towerTool = towers +createClusters.nEtaWindow = 9 +createClusters.nPhiWindow = 17 +createClusters.nEtaPosition = 5 +createClusters.nPhiPosition = 11 +createClusters.nEtaDuplicates = 7 +createClusters.nPhiDuplicates = 13 +createClusters.nEtaFinal = 9 +createClusters.nPhiFinal = 17 +# Minimal energy to create a cluster in GeV (FCC-ee detectors have to reconstruct low energy particles) +createClusters.energyThreshold = 0.1 +createClusters.attachCells = True +createClusters.clusters.Path = "CaloClusters" +createClusters.clusterCells.Path = "CaloClusterCells" +ApplicationMgr().TopAlg += [createClusters] + + +################ Output +from Configurables import PodioOutput +out = PodioOutput("out") +out.outputCommands = ["keep *", "drop ECalBarrelHits", "drop HCal*", "drop ECalBarrelCellsStep*", "drop ECalBarrelPositionedHits", "drop ECalEndcapPositionedHits", "drop emptyCaloCells", "drop CaloClusterCells"] +import uuid +out.filename = "output_fullCalo_SimAndDigi.root" +ApplicationMgr().TopAlg += [out] + + +from Configurables import CreateCaloJet + +createJets = CreateCaloJet("createJets", + InputCollection=createClusters.clusters.Path, + OutputCollection = "Jets") +# JetAlg = "antikt", +# JetRadius = 0.4, +# OutputLevel=INFO +# ) + + +ApplicationMgr().TopAlg += [createJets] + + +#CPU information +from Configurables import AuditorSvc, ChronoAuditor +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +gen.AuditExecute = True +hepmc_converter.AuditExecute = True +geantsim.AuditExecute = True +createEcalBarrelCellsStep1.AuditExecute = True +resegmentEcalBarrel.AuditExecute = True +createEcalBarrelCells.AuditExecute = True +#createHcalBarrelCells.AuditExecute = True +out.AuditExecute = True +#ApplicationMgr().ExtSvc += [audsvc] + +from Configurables import EventCounter +event_counter = EventCounter('event_counter') +event_counter.Frequency = 10 +ApplicationMgr().TopAlg += [event_counter] + diff --git a/RecCalorimeter/tests/options/runTruthJetReco.py b/RecCalorimeter/tests/options/runTruthJetReco.py new file mode 100644 index 00000000..d1647521 --- /dev/null +++ b/RecCalorimeter/tests/options/runTruthJetReco.py @@ -0,0 +1,119 @@ +import os +from Gaudi.Configuration import * +from GaudiKernel import SystemOfUnits as units +from GaudiKernel import PhysicalConstants as constants +from GaudiKernel.SystemOfUnits import MeV, GeV, tesla + +from Configurables import ApplicationMgr +ApplicationMgr().EvtSel = 'NONE' +ApplicationMgr().EvtMax = 2 +ApplicationMgr().OutputLevel = INFO +ApplicationMgr().StopOnSignal = True +ApplicationMgr().ExtSvc += ['RndmGenSvc'] + +from Configurables import k4DataSvc +## Data service +podioevent = k4DataSvc("EventDataSvc") +ApplicationMgr().ExtSvc += [podioevent] + +from Configurables import MomentumRangeParticleGun +guntool = MomentumRangeParticleGun() +guntool.ThetaMin = 45 * constants.pi / 180. +guntool.ThetaMax = 135 * constants.pi / 180. +guntool.PhiMin = 0. +guntool.PhiMax = 2. * constants.pi +guntool.MomentumMin = 1. *units.GeV +guntool.MomentumMax = 1. *units.GeV +guntool.PdgCodes = [11] # 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+ + +from Configurables import GenAlg +gen = GenAlg() +gen.SignalProvider=guntool +gen.hepmc.Path = "hepmc" +ApplicationMgr().TopAlg += [gen] + +from Configurables import HepMCToEDMConverter +## reads an HepMC::GenEvent from the data service and writes a collection of EDM Particles +hepmc_converter = HepMCToEDMConverter("Converter") +hepmc_converter.hepmc.Path="hepmc" +hepmc_converter.GenParticles.Path="GenParticles" +ApplicationMgr().TopAlg += [hepmc_converter] + + + +################## Simulation setup +# Detector geometry +from Configurables import GeoSvc +geoservice = GeoSvc("GeoSvc") +# if FCC_DETECTORS is empty, this should use relative path to working directory +path_to_detectors = os.environ.get("K4GEO", "") +detectors = [ + 'FCCee/ALLEGRO/compact/ALLEGRO_o1_v01/ALLEGRO_o1_v01.xml' +] +# prefix all xmls with path_to_detectors +for det in detectors: + geoservice.detectors += [os.path.join(path_to_detectors, det)] +geoservice.OutputLevel = INFO +ApplicationMgr().ExtSvc += [geoservice] + + +# Geant4 service +# Configures the Geant simulation: geometry, physics list and user actions +from Configurables import SimG4Svc +# giving the names of tools will initialize the tools of that type +geantservice = SimG4Svc("SimG4Svc") +geantservice.detector = "SimG4DD4hepDetector" +geantservice.physicslist = "SimG4FtfpBert" +geantservice.actions = "SimG4FullSimActions" +# Fixed seed to have reproducible results, change it for each job if you split one production into several jobs +# Mind that if you leave Gaudi handle random seed and some job start within the same second (very likely) you will have duplicates +geantservice.randomNumbersFromGaudi = False +geantservice.seedValue = 4242 +# Range cut +geantservice.g4PreInitCommands += ["/run/setCut 0.1 mm"] +ApplicationMgr().ExtSvc += [geantservice] + + + +# Magnetic field +from Configurables import SimG4ConstantMagneticFieldTool +field = SimG4ConstantMagneticFieldTool("SimG4ConstantMagneticFieldTool") +field.FieldComponentZ = -2 * units.tesla +field.FieldOn = True +field.IntegratorStepper="ClassicalRK4" + +# Geant4 algorithm +# Translates EDM to G4Event, passes the event to G4, writes out outputs via tools +# and a tool that saves the calorimeter hits + + +from Configurables import CreateTruthJet +from Configurables import FilterTruthParticlesForGenJets + +filterMCParticles = FilterTruthParticlesForGenJets("filterMCParticles", + InputCollection="GenParticles", + OutputCollection = "FilteredGenParticles") +ApplicationMgr().TopAlg += [filterMCParticles] + +createJets = CreateTruthJet("createTruthJets", + InputCollection="FilteredGenParticles", + OutputCollection = "TruthJets") + + +ApplicationMgr().TopAlg += [createJets] + + +#CPU information +from Configurables import AuditorSvc, ChronoAuditor +chra = ChronoAuditor() +audsvc = AuditorSvc() +audsvc.Auditors = [chra] +gen.AuditExecute = True +hepmc_converter.AuditExecute = True +#out.AuditExecute = True + +from Configurables import EventCounter +event_counter = EventCounter('event_counter') +event_counter.Frequency = 10 +ApplicationMgr().TopAlg += [event_counter] + diff --git a/cmake/FindFastJet.cmake b/cmake/FindFastJet.cmake new file mode 100644 index 00000000..d4580ab2 --- /dev/null +++ b/cmake/FindFastJet.cmake @@ -0,0 +1,29 @@ +# - Locate FastJet library +# Defines: +# +# FASTJET_FOUND +# FASTJET_INCLUDE_DIR +# FASTJET_INCLUDE_DIRS (not cached) +# FASTJET_LIBRARY +# FASTJET_LIBRARIES (not cached) +# FASTJET_LIBRARY_DIRS (not cached) + +find_path(FASTJET_INCLUDE_DIR fastjet/version.hh + HINTS $ENV{FASTJET_ROOT_DIR}/include ${FASTJET_ROOT_DIR}/include) + +find_library(FASTJET_LIBRARY NAMES fastjet + HINTS $ENV{FASTJET_ROOT_DIR}/lib ${FASTJET_ROOT_DIR}/lib) + +find_library(FASTJETPLUGINS_LIBRARY NAMES fastjetplugins + HINTS $ENV{FASTJET_ROOT_DIR}/lib ${FASTJET_ROOT_DIR}/lib) + +# handle the QUIETLY and REQUIRED arguments and set FASTJET_FOUND to TRUE if +# all listed variables are TRUE +INCLUDE(FindPackageHandleStandardArgs) +FIND_PACKAGE_HANDLE_STANDARD_ARGS(FastJet DEFAULT_MSG FASTJET_INCLUDE_DIR FASTJET_LIBRARY) + +mark_as_advanced(FASTJET_FOUND FASTJET_INCLUDE_DIR FASTJET_LIBRARY) + +set(FASTJET_INCLUDE_DIRS ${FASTJET_INCLUDE_DIR}) +set(FASTJET_LIBRARIES ${FASTJET_LIBRARY} ${FASTJETPLUGINS_LIBRARY}) +get_filename_component(FASTJET_LIBRARY_DIRS ${FASTJET_LIBRARY} PATH)