diff --git a/analyzers/dataframe/FCCAnalyses/Association.h b/analyzers/dataframe/FCCAnalyses/Association.h new file mode 100644 index 00000000000..2c28725745e --- /dev/null +++ b/analyzers/dataframe/FCCAnalyses/Association.h @@ -0,0 +1,72 @@ +#ifndef ASSOCIATION_ANALYZERS_H +#define ASSOCIATION_ANALYZERS_H + +// EDM4hep +#include "edm4hep/MCRecoParticleAssociationCollection.h" + + +namespace FCCAnalyses :: Association { + /** + * \brief Analyzer to select associations with a mcparticle with a + * specified PDG ID. + * + * \param pdg Desired PDG ID of the partile. + */ + struct sel_PDG { + sel_PDG(const int pdg); + const int m_pdg; + template + T operator() ( + const T& inAssocColl); + }; + + + /** + * \brief Analyzer to select associations with a mcparticle with a + * specified absolute value of PDG ID. + * + * \param pdg Desired absolute value of PDG ID of the partile. + */ + struct sel_absPDG { + const int m_pdg; + + sel_absPDG(const int pdg) : m_pdg(pdg) { + if (m_pdg < 0) { + throw std::invalid_argument( + "Association::sel_absPDG: Received negative value!"); + } + }; + + template + auto operator() (const T& inAssocColl) { + T result; + result.setSubsetCollection(); + + for (const auto& assoc: inAssocColl) { + const auto& particle = assoc.getSim(); + if (std::abs(particle.getPDG()) == m_pdg) { + result.push_back(assoc); + } + } + + return result; + }; + }; + + + /** + * \brief Analyzer to select associations with a mcparticle with a + * specified generator status. + * + * \param pdg Desired generator status of the partile. + */ + struct sel_genStatus { + sel_genStatus(const int status); + const int m_status; + template + T operator() ( + const T& inAssocColl); + }; +} /* FCCAnalyses::Association */ + +#endif /* ASSOCIATION_ANALYZERS_H */ diff --git a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle.h b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle.h index 9d2c1d9630a..ade8e44d1df 100644 --- a/analyzers/dataframe/FCCAnalyses/ReconstructedParticle.h +++ b/analyzers/dataframe/FCCAnalyses/ReconstructedParticle.h @@ -55,8 +55,8 @@ namespace ReconstructedParticle{ struct sel_absType { sel_absType(const int type); const int m_type; - ROOT::VecOps::RVec - operator()(ROOT::VecOps::RVec in); + edm4hep::ReconstructedParticleCollection + operator()(edm4hep::ReconstructedParticleCollection& in); }; /// select ReconstructedParticles with transverse momentum greater than a minimum value [GeV] @@ -134,6 +134,19 @@ namespace ReconstructedParticle{ }; + /** + * \brief Analyzer to select reconstructed particles by generator status + * + * \param status Desired generator status of the particles + */ + struct sel_genStatus { + sel_genStatus(const int status); + const int m_status; + edm4hep::ReconstructedParticleCollection operator() ( + const edm4hep::MCRecoParticleAssociationCollection& inAssocColl); + }; + + /// return reconstructed particles ROOT::VecOps::RVec get(ROOT::VecOps::RVec index, ROOT::VecOps::RVec in); diff --git a/analyzers/dataframe/src/Association.cc b/analyzers/dataframe/src/Association.cc new file mode 100644 index 00000000000..0e113a5d4b9 --- /dev/null +++ b/analyzers/dataframe/src/Association.cc @@ -0,0 +1,68 @@ +#include "FCCAnalyses/Association.h" + +// std +#include +#include + + +namespace FCCAnalyses :: Association { + sel_PDG::sel_PDG(const int pdg) : m_pdg(pdg) {}; + + template + T sel_PDG::operator() (const T& inAssocColl) { + T result; + result.setSubsetCollection(); + + for (const auto& assoc: inAssocColl) { + const auto& particle = assoc.getSim(); + if (particle.getPDG() == m_pdg) { + result.push_back(assoc); + } + } + + return result; + } + + + /* + sel_absPDG::sel_absPDG(const int pdg) : m_pdg(pdg) { + if (m_pdg < 0) { + throw std::invalid_argument( + "Association::sel_absPDG: Received negative value!"); + } + } + + template + auto sel_absPDG::operator() (const T& inAssocColl) { + T result; + result.setSubsetCollection(); + + for (const auto& assoc: inAssocColl) { + const auto& particle = assoc.getSim(); + if (std::abs(particle.getPDG()) == m_pdg) { + result.push_back(assoc); + } + } + + return result; + } + */ + + + sel_genStatus::sel_genStatus(const int status) : m_status(status) {}; + + template + T sel_genStatus::operator() (const T& inAssocColl) { + T result; + result.setSubsetCollection(); + + for (const auto& assoc: inAssocColl) { + const auto& particle = assoc.getSim(); + if (particle.getGeneratorStatus() == m_status) { + result.push_back(assoc); + } + } + + return result; + } +} /* FCCAnalyses::Association */ diff --git a/analyzers/dataframe/src/ReconstructedParticle.cc b/analyzers/dataframe/src/ReconstructedParticle.cc index 8e774002952..146c9570854 100644 --- a/analyzers/dataframe/src/ReconstructedParticle.cc +++ b/analyzers/dataframe/src/ReconstructedParticle.cc @@ -5,6 +5,14 @@ #include #include +// ROOT +#include +#include + + +#define rdfInfo R__LOG_INFO(ROOT::Detail::RDF::RDFLogChannel()) + + namespace FCCAnalyses{ namespace ReconstructedParticle{ @@ -24,6 +32,7 @@ ROOT::VecOps::RVec sel_type::operator()( return result; } + /// sel_absType sel_absType::sel_absType(const int type) : m_type(type) { if (m_type < 0) { @@ -32,18 +41,21 @@ sel_absType::sel_absType(const int type) : m_type(type) { } } -ROOT::VecOps::RVec sel_absType::operator()( - ROOT::VecOps::RVec in) { - ROOT::VecOps::RVec result; - result.reserve(in.size()); - for (size_t i = 0; i < in.size(); ++i) { - if (std::abs(in[i].type) == m_type) { - result.emplace_back(in[i]); +edm4hep::ReconstructedParticleCollection sel_absType::operator()( + edm4hep::ReconstructedParticleCollection& inColl) { + edm4hep::ReconstructedParticleCollection result; + result.setSubsetCollection(); + + for (const auto& particle: inColl) { + if (std::abs(particle.getType()) == m_type) { + result.push_back(particle); } } + return result; } + /// sel_pt sel_pt::sel_pt(float arg_min_pt) : m_min_pt(arg_min_pt) {}; ROOT::VecOps::RVec sel_pt::operator() (ROOT::VecOps::RVec in) { @@ -113,6 +125,8 @@ edm4hep::ReconstructedParticleCollection selPDG::operator() ( edm4hep::ReconstructedParticleCollection result; result.setSubsetCollection(); + rdfInfo << "Assoc coll size: " << inAssocColl.size(); + for (const auto& assoc: inAssocColl) { const auto& particle = assoc.getSim(); if (m_chargeConjugateAllowed) { @@ -148,6 +162,40 @@ edm4hep::ReconstructedParticleCollection selUpTo::operator() ( return result; } +sel_genStatus::sel_genStatus(int status) : m_status(status) {}; + +edm4hep::ReconstructedParticleCollection sel_genStatus::operator() ( + const edm4hep::MCRecoParticleAssociationCollection& inAssocColl) { + edm4hep::ReconstructedParticleCollection result; + result.setSubsetCollection(); + + + // std::cout << "******************************************" << std::endl; + for (const auto& assoc: inAssocColl) { + const auto& mcParticle = assoc.getSim(); + const auto& recoParticle = assoc.getRec(); + // std::cout << "------------------------------------" << std::endl; + // std::cout << "gen status: " << mcParticle.getGeneratorStatus() << std::endl; + // std::cout << "energy: " << mcParticle.getEnergy() << std::endl; + // std::cout << "pdg: " << mcParticle.getPDG() << std::endl; + // std::cout << "" << std::endl; + // std::cout << "type: " << recoParticle.getType() << std::endl; + // std::cout << "energy: " << recoParticle.getEnergy() << std::endl; + // std::cout << "momentum, x: " << recoParticle.getMomentum().x << std::endl; + // std::cout << "momentum, y: " << recoParticle.getMomentum().y << std::endl; + // std::cout << "momentum, z: " << recoParticle.getMomentum().z << std::endl; + // std::cout << "position, x: " << recoParticle.getReferencePoint().x << std::endl; + // std::cout << "position, y: " << recoParticle.getReferencePoint().y << std::endl; + // std::cout << "position, z: " << recoParticle.getReferencePoint().z << std::endl; + // std::cout << "goodnessOfPID: " << recoParticle.getGoodnessOfPID() << std::endl; + if (mcParticle.getGeneratorStatus() == m_status) { + result.push_back(assoc.getRec()); + } + } + + return result; +} + resonanceBuilder::resonanceBuilder(float arg_resonance_mass) {m_resonance_mass = arg_resonance_mass;} ROOT::VecOps::RVec resonanceBuilder::operator()(ROOT::VecOps::RVec legs) { diff --git a/e4hsource/src/EDM4hepDataSource.cxx b/e4hsource/src/EDM4hepDataSource.cxx index ba551397843..5edc8bf770b 100644 --- a/e4hsource/src/EDM4hepDataSource.cxx +++ b/e4hsource/src/EDM4hepDataSource.cxx @@ -42,7 +42,7 @@ namespace FCCAnalyses { * \brief Setup input for the EDM4hepDataSource. */ void EDM4hepDataSource::SetupInput(int nEvents) { - std::cout << "EDM4hepDataSource: Constructing the source ..." << std::endl; + // std::cout << "EDM4hepDataSource: Constructing the source ..." << std::endl; if (m_filePathList.empty()) { throw std::runtime_error("EDM4hepDataSource: No input files provided!"); @@ -75,11 +75,13 @@ namespace FCCAnalyses { // Determine over how many events to run if (nEventsInFiles > 0) { + /* std::cout << "EDM4hepDataSource: Found " << nEventsInFiles << " events in files: \n"; for (const auto& filePath : m_filePathList) { std::cout << " - " << filePath << "\n"; } + */ } else { throw std::runtime_error("EDM4hepDataSource: No events found!"); } @@ -96,18 +98,18 @@ namespace FCCAnalyses { m_nEvents = nEventsInFiles; } - std::cout << "EDM4hepDataSource: Running over " << m_nEvents << " events." - << std::endl; + // std::cout << "EDM4hepDataSource: Running over " << m_nEvents << " events." + // << std::endl; // Get collections stored in the files std::vector collNames = frame.getAvailableCollections(); - std::cout << "EDM4hepDataSource: Found following collections:\n"; + // std::cout << "EDM4hepDataSource: Found following collections:\n"; for (auto& collName: collNames) { const podio::CollectionBase* coll = frame.get(collName); if (coll->isValid()) { m_columnNames.emplace_back(collName); m_columnTypes.emplace_back(coll->getValueTypeName()); - std::cout << " - " << collName << "\n"; + // std::cout << " - " << collName << "\n"; } } } @@ -118,8 +120,8 @@ namespace FCCAnalyses { */ void EDM4hepDataSource::SetNSlots(unsigned int nSlots) { - std::cout << "EDM4hepDataSource: Setting num. of slots to: " << nSlots - << std::endl; + // std::cout << "EDM4hepDataSource: Setting num. of slots to: " << nSlots + // << std::endl; m_nSlots = nSlots; if (m_nSlots > m_nEvents) { @@ -160,8 +162,7 @@ namespace FCCAnalyses { */ void EDM4hepDataSource::Initialize() { - std::cout << "EDM4hepDataSource: Initializing the source ..." << std::endl; - + // std::cout << "EDM4hepDataSource: Initializing the source ..." << std::endl; } @@ -171,7 +172,7 @@ namespace FCCAnalyses { */ std::vector> EDM4hepDataSource::GetEntryRanges() { - std::cout << "EDM4hepDataSource: Getting entry ranges ..." << std::endl; + // std::cout << "EDM4hepDataSource: Getting entry ranges ..." << std::endl; std::vector> rangesToBeProcessed; for (auto& range: m_rangesAvailable) { @@ -192,6 +193,7 @@ namespace FCCAnalyses { } + /* std::cout << "EDM4hepDataSource: Ranges to be processed:\n"; for (auto& range: rangesToBeProcessed) { std::cout << " {" << range.first << ", " << range.second @@ -208,6 +210,7 @@ namespace FCCAnalyses { } else { std::cout << "EDM4hepDataSource: No more remaining ranges.\n"; } + */ return rangesToBeProcessed; } @@ -219,8 +222,8 @@ namespace FCCAnalyses { */ void EDM4hepDataSource::InitSlot(unsigned int slot, ULong64_t firstEntry) { - std::cout << "EDM4hepDataSource: Initializing slot: " << slot - << " with first entry " << firstEntry << std::endl; + // std::cout << "EDM4hepDataSource: Initializing slot: " << slot + // << " with first entry " << firstEntry << std::endl; } @@ -230,8 +233,8 @@ namespace FCCAnalyses { */ bool EDM4hepDataSource::SetEntry(unsigned int slot, ULong64_t entry) { - std::cout << "EDM4hepDataSource: In slot: " << slot << ", setting entry: " - << entry << std::endl; + // std::cout << "EDM4hepDataSource: In slot: " << slot << ", setting entry: " + // << entry << std::endl; m_frames[slot] = std::make_unique( podio::Frame(m_podioReaders[slot]->readEntry("events", entry))); @@ -239,12 +242,14 @@ namespace FCCAnalyses { for (auto& collectionIndex: m_activeCollections) { m_Collections[collectionIndex][slot] = m_frames[slot]->get(m_columnNames.at(collectionIndex)); + /* std::cout << "CollName: " << m_columnNames.at(collectionIndex) << "\n"; std::cout << "Address: " << m_Collections[collectionIndex][slot] << "\n"; std::cout << "Coll size: " << m_Collections[collectionIndex][slot]->size() << "\n"; if (m_Collections[collectionIndex][slot]->isValid()) { std::cout << "Collection valid\n"; } + */ } return true; @@ -257,7 +262,7 @@ namespace FCCAnalyses { */ void EDM4hepDataSource::FinalizeSlot(unsigned int slot) { - std::cout << "EDM4hepDataSource: Finalizing slot: " << slot << std::endl; + // std::cout << "EDM4hepDataSource: Finalizing slot: " << slot << std::endl; // std::cout << "Reader: " << &m_podioReaderRefs[slot].get() << std::endl; // for (auto& collectionIndex: m_activeCollections) { @@ -276,7 +281,7 @@ namespace FCCAnalyses { */ void EDM4hepDataSource::Finalize() { - std::cout << "EDM4hepDataSource: Finalizing ..." << std::endl; + // std::cout << "EDM4hepDataSource: Finalizing ..." << std::endl; } @@ -287,9 +292,11 @@ namespace FCCAnalyses { Record_t EDM4hepDataSource::GetColumnReadersImpl(std::string_view columnName, const std::type_info& typeInfo) { + /* std::cout << "EDM4hepDataSource: Getting column reader implementation for column:\n" << " " << columnName << "\n with type: " << typeInfo.name() << std::endl; + */ auto itr = std::find(m_columnNames.begin(), m_columnNames.end(), columnName); @@ -301,12 +308,14 @@ namespace FCCAnalyses { } auto columnIndex = std::distance(m_columnNames.begin(), itr); m_activeCollections.emplace_back(columnIndex); + /* std::cout << "EDM4hepDataSource: Active collections so far:\n" << " "; for (auto& i: m_activeCollections) { std::cout << i << ", "; } std::cout << std::endl; + */ Record_t columnReaders(m_nSlots); for (size_t slotIndex = 0; slotIndex < m_nSlots; ++slotIndex) { @@ -327,7 +336,7 @@ namespace FCCAnalyses { */ const std::vector& EDM4hepDataSource::GetColumnNames() const { - std::cout << "EDM4hepDataSource: Looking for column names" << std::endl; + // std::cout << "EDM4hepDataSource: Looking for column names" << std::endl; return m_columnNames; } @@ -337,8 +346,8 @@ namespace FCCAnalyses { */ bool EDM4hepDataSource::HasColumn(std::string_view columnName) const { - std::cout << "EDM4hepDataSource: Looking for column: " << columnName - << std::endl; + // std::cout << "EDM4hepDataSource: Looking for column: " << columnName + // << std::endl; if (std::find(m_columnNames.begin(), m_columnNames.end(), @@ -355,15 +364,15 @@ namespace FCCAnalyses { */ std::string EDM4hepDataSource::GetTypeName(std::string_view columnName) const { - std::cout << "EDM4hepDataSource: Looking for type name of column: " - << columnName << std::endl; + // std::cout << "EDM4hepDataSource: Looking for type name of column: " + // << columnName << std::endl; auto itr = std::find(m_columnNames.begin(), m_columnNames.end(), columnName); if (itr != m_columnNames.end()) { auto i = std::distance(m_columnNames.begin(), itr); - std::cout << "EDM4hepDataSource: Found type name: " - << m_columnTypes.at(i) << std::endl; + // std::cout << "EDM4hepDataSource: Found type name: " + // << m_columnTypes.at(i) << std::endl; return m_columnTypes.at(i) + "Collection"; } diff --git a/tests/unittest/ReconstructedParticle.cpp b/tests/unittest/ReconstructedParticle.cpp index d68e19de480..c5bf47445c4 100644 --- a/tests/unittest/ReconstructedParticle.cpp +++ b/tests/unittest/ReconstructedParticle.cpp @@ -25,24 +25,24 @@ TEST_CASE("sel_type", "[ReconstructedParticle]") { } TEST_CASE("sel_absType", "[ReconstructedParticle]") { - ROOT::VecOps::RVec pVec; - edm4hep::ReconstructedParticleData p1; - p1.type = 11; - pVec.push_back(p1); - edm4hep::ReconstructedParticleData p2; - p2.type = 13; - pVec.push_back(p2); - edm4hep::ReconstructedParticleData p3; - p3.type = -11; - pVec.push_back(p3); - edm4hep::ReconstructedParticleData p4; - p4.type = -13; - pVec.push_back(p4); + edm4hep::ReconstructedParticleCollection pColl; + edm4hep::MutableReconstructedParticle p1; + p1.setType(11); + pColl.push_back(p1); + edm4hep::MutableReconstructedParticle p2; + p2.setType(13); + pColl.push_back(p2); + edm4hep::MutableReconstructedParticle p3; + p3.setType(-11); + pColl.push_back(p3); + edm4hep::MutableReconstructedParticle p4; + p4.setType(-13); + pColl.push_back(p4); FCCAnalyses::ReconstructedParticle::sel_absType selAbsType{11}; - auto res = selAbsType(pVec); + auto res = selAbsType(pColl); REQUIRE(res.size() == 2); - REQUIRE(res[0].type == 11); - REQUIRE(res[1].type == -11); + REQUIRE(res[0].getType() == 11); + REQUIRE(res[1].getType() == -11); } TEST_CASE("sel_absType__neg_type", "[ReconstructedParticle]") {