Jets/src/TruthJetAlgorithm.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Jets/TruthJetAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Utilities/Logger.hpp"
0013 #include "ActsExamples/EventData/SimParticle.hpp"
0014 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0015 #include "ActsFatras/EventData/ProcessType.hpp"
0016
0017 #include <ostream>
0018 #include <stdexcept>
0019
0020 #include <fastjet/ClusterSequence.hh>
0021 #include <fastjet/JetDefinition.hh>
0022 #include <fastjet/PseudoJet.hh>
0023
0024 namespace ActsExamples {
0025
0026 TruthJetAlgorithm::TruthJetAlgorithm(const Config& cfg,
0027                                      Acts::Logging::Level lvl)
0028     : IAlgorithm("TruthJetAlgorithm", lvl), m_cfg(cfg) {
0029   if (m_cfg.inputTruthParticles.empty()) {
0030     throw std::invalid_argument("Input particles collection is not configured");
0031   }
0032   m_inputTruthParticles.initialize(m_cfg.inputTruthParticles);
0033   m_outputJets.initialize(m_cfg.outputJets);
0034 }
0035
0036 ProcessCode ActsExamples::TruthJetAlgorithm::execute(
0037     const ActsExamples::AlgorithmContext& ctx) const {
0038   const auto& truthParticles = m_inputTruthParticles(ctx);
0039
0040   ACTS_DEBUG("Number of truth particles: " << truthParticles.size());
0041
0042   const fastjet::JetDefinition defaultJetDefinition =
0043       fastjet::JetDefinition(fastjet::antikt_algorithm, 0.4);
0044
0045   // Get the 4-momentum information from the simulated truth particles
0046   // and create fastjet::PseudoJet objects
0047   std::vector<fastjet::PseudoJet> inputPseudoJets;
0048
0049   int particleIndex = 0;
0050   for (const auto& particle : truthParticles) {
0051     fastjet::PseudoJet pseudoJet(particle.momentum().x(),
0052                                  particle.momentum().y(),
0053                                  particle.momentum().z(), particle.energy());
0054
0055     pseudoJet.set_user_index(particleIndex);
0056     inputPseudoJets.push_back(pseudoJet);
0057     particleIndex++;
0058   }
0059   ACTS_DEBUG("Number of input pseudo jets: " << inputPseudoJets.size());
0060   // Run the jet clustering
0061   fastjet::ClusterSequence clusterSeq(inputPseudoJets, defaultJetDefinition);
0062   // Get the jets above a certain pt threshold
0063   std::vector<fastjet::PseudoJet> jets =
0064       sorted_by_pt(clusterSeq.inclusive_jets(m_cfg.jetPtMin));
0065   ACTS_DEBUG("Number of clustered jets: " << jets.size());
0066   // Store the jets in the output data handle
0067   m_outputJets(ctx, std::move(jets));
0068
0069   return ProcessCode::SUCCESS;
0070 }
0071
0072 ProcessCode ActsExamples::TruthJetAlgorithm::finalize() {
0073   ACTS_INFO("Finalizing truth jet clustering");
0074   return ProcessCode::SUCCESS;
0075 }
0076
0077 };  // namespace ActsExamples