diff --git a/include/ral/LogicalOperators.h b/include/ral/LogicalOperators.h index 7466608..d2fb890 100644 --- a/include/ral/LogicalOperators.h +++ b/include/ral/LogicalOperators.h @@ -1,8 +1,15 @@ #ifndef LOGICAL_OPERATORS_RAL_H #define LOGICAL_OPERATORS_RAL_H +// std +#include + +// ROOT #include "ROOT/RVec.hxx" +// PODIO +#include "podio/CollectionBase.h" + namespace k4::ral { namespace LogicalOperators { @@ -10,8 +17,43 @@ namespace LogicalOperators { enum class ComparisonOperator { LESS, LESSEQ, EQ, GREATEREQ, GREATER }; template -ROOT::VecOps::RVec filter(ROOT::VecOps::RVec mask, - ROOT::VecOps::RVec collection); +ROOT::VecOps::RVec filter(const ROOT::VecOps::RVec &mask, + const ROOT::VecOps::RVec &collection) { + if (mask.size() != collection.size()) { + auto msg = "Different vector lengths: " + std::to_string(mask.size()) + + " vs. " + std::to_string(collection.size()) + "!"; + throw std::length_error(msg); + } + + ROOT::VecOps::RVec result; + for (int i = 0; i < collection.size(); i++) { + if (mask[i]) { + result.emplace_back(collection[i]); + } + } + + return result; +} + +template ::value>> +C filter(const ROOT::VecOps::RVec &mask, const C &collection) { + if (mask.size() != collection.size()) { + auto msg = "Different vector lengths: " + std::to_string(mask.size()) + + " vs. " + std::to_string(collection.size()) + "!"; + throw std::length_error(msg); + } + + C result; + result.setSubsetCollection(); + for (int i = 0; i < collection.size(); i++) { + if (mask[i]) { + result.push_back(collection[i]); + } + } + + return result; +} ROOT::VecOps::RVec operator&&(const ROOT::VecOps::RVec &vec1, const ROOT::VecOps::RVec &vec2); diff --git a/include/ral/ReconstructedParticle.h b/include/ral/ReconstructedParticle.h index c3bfec1..da744f6 100644 --- a/include/ral/ReconstructedParticle.h +++ b/include/ral/ReconstructedParticle.h @@ -5,7 +5,7 @@ #include "Math/Vector4D.h" #include "ROOT/RVec.hxx" #include "edm4hep/ClusterData.h" -#include "edm4hep/ReconstructedParticleData.h" +#include "edm4hep/ReconstructedParticleCollection.h" #include "edm4hep/TrackData.h" #include "ral/LogicalOperators.h" #include @@ -218,6 +218,15 @@ get_absq(ROOT::VecOps::RVec particles); ROOT::VecOps::RVec get_m(ROOT::VecOps::RVec particles); +/** + * Get mass member from ReconstructedParticles + * + * @param particles List of reconstructed particles in an event + * + */ +ROOT::VecOps::RVec +get_m(const edm4hep::ReconstructedParticleCollection &particles); + /** * Get goodnessOfPID member from ReconstructedParticles * @@ -586,6 +595,18 @@ ROOT::VecOps::RVec sel_e(LogicalOperators::ComparisonOperator op, float value, ROOT::VecOps::RVec particles); +/** + * Select a subcollection of ReconstructedParticles based on the energy + * + * @param op Comparison operator to apply with the value + * @param value Value of the energy that is used in the comparison + * @param particles Collection of reconstructed particles in an event + * + */ +edm4hep::ReconstructedParticleCollection +sel_e(LogicalOperators::ComparisonOperator op, float value, + const edm4hep::ReconstructedParticleCollection &particles); + /** * Select a subcollection of ReconstructedParticles based on the 3D momentum * modulus diff --git a/src/LogicalOperators.cc b/src/LogicalOperators.cc index 8838fd8..a2ac19c 100644 --- a/src/LogicalOperators.cc +++ b/src/LogicalOperators.cc @@ -1,6 +1,5 @@ #include "ral/LogicalOperators.h" #include "ROOT/RVec.hxx" -#include #include #include @@ -8,30 +7,6 @@ namespace k4::ral { namespace LogicalOperators { -template -ROOT::VecOps::RVec filter(ROOT::VecOps::RVec mask, - ROOT::VecOps::RVec collection) { - if (mask.size() != collection.size()) { - auto msg = "Different vector lenghts: " + std::to_string(mask.size()) + - " and " + std::to_string(collection.size()) + "."; - throw std::length_error(msg); - } - ROOT::VecOps::RVec result; - for (int i = 0; i < collection.size(); i++) { - if (mask[i]) { - result.emplace_back(collection[i]); - } - } - return result; -} - -/// It is important to instanciate the methods that we want to use in order -/// them to be available -template ROOT::VecOps::RVec -filter( - ROOT::VecOps::RVec mask, - ROOT::VecOps::RVec collection); - ROOT::VecOps::RVec operator&&(const ROOT::VecOps::RVec &vec1, const ROOT::VecOps::RVec &vec2) { if (vec1.size() != vec2.size()) { diff --git a/src/ReconstructedParticle.cc b/src/ReconstructedParticle.cc index 4b14195..221e287 100644 --- a/src/ReconstructedParticle.cc +++ b/src/ReconstructedParticle.cc @@ -49,6 +49,16 @@ get_e(ROOT::VecOps::RVec particles) { return result; } +ROOT::VecOps::RVec +get_e(const edm4hep::ReconstructedParticleCollection &particles) { + ROOT::VecOps::RVec result; + result.reserve(particles.size()); + for (edm4hep::ReconstructedParticle p : particles) { + result.emplace_back(p.getEnergy()); + } + return result; +} + ROOT::VecOps::RVec get_p(ROOT::VecOps::RVec particles) { ROOT::VecOps::RVec result; @@ -247,6 +257,16 @@ get_m(ROOT::VecOps::RVec particles) { return result; } +ROOT::VecOps::RVec +get_m(const edm4hep::ReconstructedParticleCollection &particles) { + ROOT::VecOps::RVec result; + result.reserve(particles.size()); + for (const auto &p : particles) { + result.emplace_back(p.getMass()); + } + return result; +} + ROOT::VecOps::RVec get_goodnessOfPID( ROOT::VecOps::RVec particles) { ROOT::VecOps::RVec result; @@ -420,7 +440,7 @@ int print_goodnessOfPID::operator()( #define MASKING(T, getter, op, value, collection, result) \ ROOT::VecOps::RVec vec = getter(collection); \ - for (T & item : vec) { \ + for (const T &item : vec) { \ switch (op) { \ case LogicalOperators::ComparisonOperator::LESS: \ result.emplace_back(item < value); \ @@ -449,6 +469,15 @@ mask_e(LogicalOperators::ComparisonOperator op, float value, return result; } +ROOT::VecOps::RVec +mask_e(LogicalOperators::ComparisonOperator op, float value, + const edm4hep::ReconstructedParticleCollection &particles) { + ROOT::VecOps::RVec result; + result.reserve(particles.size()); + MASKING(float, get_e, op, value, particles, result); + return result; +} + ROOT::VecOps::RVec mask_pmod(LogicalOperators::ComparisonOperator op, float value, ROOT::VecOps::RVec particles) { @@ -631,6 +660,13 @@ sel_e(LogicalOperators::ComparisonOperator op, float value, mask, particles); } +edm4hep::ReconstructedParticleCollection +sel_e(LogicalOperators::ComparisonOperator op, float value, + const edm4hep::ReconstructedParticleCollection &particles) { + auto mask = mask_e(op, value, particles); + return LogicalOperators::filter(mask, particles); +} + ROOT::VecOps::RVec sel_pmod(LogicalOperators::ComparisonOperator op, float value, ROOT::VecOps::RVec particles) { diff --git a/test/integration/CMakeLists.txt b/test/integration/CMakeLists.txt index cb05153..87f6053 100644 --- a/test/integration/CMakeLists.txt +++ b/test/integration/CMakeLists.txt @@ -1,2 +1,5 @@ add_test(NAME integration_ReconstructedParticle COMMAND python "${CMAKE_CURRENT_LIST_DIR}/ReconstructedParticle.py") + + add_test(NAME integration_ReconstructedParticleSource + COMMAND python "${CMAKE_CURRENT_LIST_DIR}/ReconstructedParticleSource.py") diff --git a/test/integration/ReconstructedParticle.py b/test/integration/ReconstructedParticle.py index 13ce386..2116acd 100644 --- a/test/integration/ReconstructedParticle.py +++ b/test/integration/ReconstructedParticle.py @@ -66,7 +66,7 @@ def get_file(url: str, path: str) -> None: if not ROOT.gSystem.Load("libral") and ROOT.gSystem.Load("libedm4hep"): print('RAL Found!') print('EDM4HEP Found!') - if ROOT.loadRal() and ROOT.edm4hep.ReconstructedParticleData(): + if ROOT.loadRal and ROOT.edm4hep.ReconstructedParticleData: print('RAL Loaded!') print('EDM4HEP Loaded!') diff --git a/test/integration/ReconstructedParticleSource.py b/test/integration/ReconstructedParticleSource.py new file mode 100644 index 0000000..1233903 --- /dev/null +++ b/test/integration/ReconstructedParticleSource.py @@ -0,0 +1,96 @@ +import sys +import subprocess +import urllib.request +from pathlib import Path +import ROOT + +def get_file(url: str, path: str) -> None: + size = 0 + with urllib.request.urlopen(url) as stream, open(path, "wb") as file: + buf_size = 1024 * 8 + while block := stream.read(buf_size): + file.write(block) + size += buf_size + if size == 0: + raise RuntimeError(f"Cannot download from url: {url}") + + +PYTHIA_CARD = { + "url": "https://raw.githubusercontent.com/HEP-FCC/FCC-config/winter2023/FCCee/Generator/Pythia8/p8_ee_WW_ecm240.cmd", + "path": "/tmp/test_pythia_card.cmd" +} +DETECTOR_CARD = { + "url": "https://raw.githubusercontent.com/HEP-FCC/FCC-config/winter2023/FCCee/Delphes/card_IDEA.tcl", + "path": "/tmp/test_detector_card.tcl" +} +OUTPUT_CARD = { + "url": "https://raw.githubusercontent.com/HEP-FCC/FCC-config/winter2023/FCCee/Delphes/card_IDEA.tcl", + "path": "/tmp/test_output_card.tcl" +} +ROOT_FILE = "/tmp/test.root" +FINAL_ROOT_FILE = "/tmp/test_final.root" + +if __name__ == "__main__": + + print("ReconstructedParticle Integration Test: PRE-TEST".center(80, "-")) + + if not Path(ROOT_FILE).exists(): + print("Downloading Pythia files") + get_file(**PYTHIA_CARD) + get_file(**DETECTOR_CARD) + get_file(**OUTPUT_CARD) + print("Running Pythia to generate root file") + subprocess.run( + [ + "DelphesPythia8_EDM4HEP", + DETECTOR_CARD["path"], + OUTPUT_CARD["path"], + PYTHIA_CARD["path"], + ROOT_FILE + ], + stdout=subprocess.DEVNULL, + stderr=subprocess.DEVNULL + ).check_returncode() + # else: + # subprocess.run( + # [ + # "podio-dump", + # ROOT_FILE + # ], + # stdout=subprocess.DEVNULL, + # stderr=subprocess.DEVNULL + # ).check_returncode() + + print("ReconstructedParticle Integration Test: START".center(80, "-")) + + if ROOT.podio.DataSource: + print('PODIO::DataSource Loaded!') + if ROOT.edm4hep.ReconstructedParticle: + print('EDM4HEP Loaded!') + if ROOT.gSystem.Load("libral") < 0: + print('Loading of the Key4hep RAL library failed!\nAborting...') + sys.exit(1) + if ROOT.loadRal: + print('Key4hep RAL Loaded!') + + print("Loading ROOT dataframe and test ral functions") + df = ROOT.podio.CreateDataFrame(ROOT_FILE) + ROOT.gInterpreter.ProcessLine("using namespace k4::ral;") + ROOT.gInterpreter.ProcessLine("using namespace LogicalOperators;") + df = ( + df + .Define("e1", + "ReconstructedParticle::sel_e(ComparisonOperator::LESS, 1., ReconstructedParticles)") + .Define("e2", + "ReconstructedParticle::sel_e(ComparisonOperator::GREATER, 5., e1)") + .Define("mass", + "ReconstructedParticle::get_m(e2)") + .Define("total_mass", + "ROOT::VecOps::Sum(mass)") + ) + print("Output test result in a new dataframe") + df.Snapshot("events", FINAL_ROOT_FILE, + [ + "mass", + "total_mass" + ]) diff --git a/test/unittest/ReconstructedParticle.cpp b/test/unittest/ReconstructedParticle.cpp index ea1e255..ce2e678 100644 --- a/test/unittest/ReconstructedParticle.cpp +++ b/test/unittest/ReconstructedParticle.cpp @@ -11,7 +11,6 @@ #include "ral/LogicalOperators.h" #include "ral/ReconstructedParticle.h" #include -#include using namespace k4::ral;