Commit ad753180 authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

Mutable datatypes in podio@0.14.1 and later

parent 6942bc8e
......@@ -23,9 +23,9 @@ if(ENABLE_CLANG_TIDY)
endif()
find_package(EICD REQUIRED)
find_package(EDM4HEP REQUIRED)
find_package(EDM4HEP 0.4.1 REQUIRED)
find_package(podio 0.11...0.14 REQUIRED)
find_package(podio 0.14.1 REQUIRED)
add_definitions("-Dpodio_VERSION_MAJOR=${podio_VERSION_MAJOR}")
add_definitions("-Dpodio_VERSION_MINOR=${podio_VERSION_MINOR}")
add_definitions("-Dpodio_VERSION_PATCH=${podio_VERSION_PATCH}")
......
......@@ -27,33 +27,31 @@ StatusCode PodioOutput::initialize() {
m_datatree = m_podioDataSvc->eventDataTree();
m_datatree->SetDirectory(m_file.get());
m_metadatatree = new TTree("metadata", "Metadata tree");
m_runMDtree = new TTree("run_metadata", "Run metadata tree");
m_evtMDtree = new TTree("evt_metadata", "Event metadata tree");
m_colMDtree = new TTree("col_metadata", "Collection metadata tree");
m_evtMDtree->Branch("evtMD", "GenericParameters", m_podioDataSvc->getProvider().eventMetaDataPtr() ) ;
m_switch = KeepDropSwitch(m_outputCommands);
return StatusCode::SUCCESS;
}
void PodioOutput::resetBranches(const std::vector<std::pair<std::string, podio::CollectionBase*>>& collections) {
for (auto& [collName, collBuffers] : collections) {
#if podio_VERSION_MAJOR == 0 && (podio_VERSION_MINOR < 13 || (podio_VERSION_MINOR == 13 && podio_VERSION_PATCH < 2))
auto data = collBuffers->getBufferAddress();
auto references = collBuffers->referenceCollections();
auto vecmembers = collBuffers->vectorMembers();
#else
auto buffers = collBuffers->getBuffers();
auto data = buffers.data;
auto references = buffers.references;
auto vecmembers = buffers.vectorMembers;
#endif
if (m_switch.isOn(collName)) {
// Reconnect branches and collections
m_datatree->SetBranchAddress(collName.c_str(), data);
auto colls = references;
if (colls) {
int j = 0;
for (auto& c : (*colls)) {
for (size_t j = 0; j < colls->size(); ++j) {
const auto brName = podio::root_utils::refBranch(collName, j);
m_datatree->SetBranchAddress(brName.c_str(), &c);
++j;
auto l_branch = m_datatree->GetBranch(brName.c_str());
l_branch->SetAddress(&(*colls)[j]);
}
}
auto colls_v = vecmembers;
......@@ -71,17 +69,15 @@ void PodioOutput::resetBranches(const std::vector<std::pair<std::string, podio::
}
void PodioOutput::createBranches(const std::vector<std::pair<std::string, podio::CollectionBase*>>& collections) {
// collectionID, collection type, subset collection
std::vector<std::tuple<int, std::string, bool>>* collectionInfo = new std::vector<std::tuple<int, std::string, bool>>();
collectionInfo->reserve(collections.size());
for (auto& [collName, collBuffers] : collections) {
#if podio_VERSION_MAJOR == 0 && (podio_VERSION_MINOR < 13 || (podio_VERSION_MINOR == 13 && podio_VERSION_PATCH < 2))
auto data = collBuffers->getBufferAddress();
auto references = collBuffers->referenceCollections();
auto vecmembers = collBuffers->vectorMembers();
#else
auto buffers = collBuffers->getBuffers();
auto data = buffers.data;
auto references = buffers.references;
auto vecmembers = buffers.vectorMembers;
#endif
const std::string className = collBuffers->getValueTypeName();
const std::string collClassName = "vector<" + className + "Data>";
......@@ -95,7 +91,7 @@ void PodioOutput::createBranches(const std::vector<std::pair<std::string, podio:
int j = 0;
for (auto& c : (*refColls)) {
const auto brName = podio::root_utils::refBranch(collName, j);
m_datatree->Branch(brName.c_str(), c);
m_datatree->Branch(brName.c_str(), c.get());
++j;
}
}
......@@ -111,10 +107,16 @@ void PodioOutput::createBranches(const std::vector<std::pair<std::string, podio:
}
}
const auto collID = m_podioDataSvc->getCollectionIDs()->collectionID(collName);
const auto collType = collBuffers->getValueTypeName() + "Collection";
collectionInfo->emplace_back(collID, std::move(collType), collBuffers->isSubsetCollection());
debug() << isOn << " Registering collection " << collClassName << " " << collName.c_str() << " containing type "
<< className << endmsg;
collBuffers->prepareForWrite();
}
m_metadatatree->Branch("CollectionTypeInfo", collectionInfo);
}
StatusCode PodioOutput::execute() {
......@@ -132,6 +134,7 @@ StatusCode PodioOutput::execute() {
debug() << "Filling DataTree .." << endmsg;
}
m_datatree->Fill();
m_evtMDtree->Fill();
return StatusCode::SUCCESS;
}
......@@ -183,6 +186,10 @@ StatusCode PodioOutput::finalize() {
m_metadatatree->Branch("gaudiConfigOptions", &config_data);
m_metadatatree->Branch("CollectionIDs", m_podioDataSvc->getCollectionIDs());
m_metadatatree->Fill();
m_colMDtree->Branch("colMD", "std::map<int,podio::GenericParameters>", m_podioDataSvc->getProvider().getColMetaDataMap() ) ;
m_colMDtree->Fill();
m_runMDtree->Branch("runMD", "std::map<int,podio::GenericParameters>", m_podioDataSvc->getProvider().getRunMetaDataMap() ) ;
m_runMDtree->Fill();
m_datatree->Write();
m_file->Write();
m_file->Close();
......
......@@ -49,6 +49,9 @@ private:
TTree* m_datatree;
/// The tree to be filled with meta data
TTree* m_metadatatree;
TTree* m_runMDtree;
TTree* m_evtMDtree;
TTree* m_colMDtree;
/// The stored collections
std::vector<podio::CollectionBase*> m_storedCollections;
};
......
......@@ -29,16 +29,10 @@ inline std::string vecBranch(const std::string& name, size_t index) {
inline void setCollectionAddresses(podio::CollectionBase* collection, const CollectionBranches& branches) {
#if podio_VERSION_MAJOR == 0 && (podio_VERSION_MINOR < 13 || (podio_VERSION_MINOR == 13 && podio_VERSION_PATCH < 2))
auto data = collection->getBufferAddress();
auto references = collection->referenceCollections();
auto vecmembers = collection->vectorMembers();
#else
auto buffers = collection->getBuffers();
auto data = buffers.data;
auto references = buffers.references;
auto vecmembers = buffers.vectorMembers;
#endif
if (data) {
branches.data->SetAddress(data);
......
......@@ -199,7 +199,7 @@ namespace Jug::Digi {
auto rawhits = m_outputHitCollection.createAndPut();
// find the hits that belong to the same group (for merging)
std::unordered_map<long long, std::vector<edm4hep::ConstSimCalorimeterHit>> merge_map;
std::unordered_map<long long, std::vector<edm4hep::SimCalorimeterHit>> merge_map;
for (const auto &ahit : *simhits) {
int64_t hid = (ahit.getCellID() & id_mask) | ref_mask;
auto it = merge_map.find(hid);
......
......@@ -57,7 +57,7 @@ public:
const auto pgen = std::hypot(pvec.x, pvec.y, pvec.z);
const auto momentum = pgen * m_gaussDist();
// make sure we keep energy consistent
using MomType = decltype(eicd::ReconstructedParticle().momentum().x);
using MomType = decltype(eicd::ReconstructedParticle().getMomentum().x);
const MomType energy = std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum);
const MomType px = p.getMomentum().x * momentum / pgen;
const MomType py = p.getMomentum().y * momentum / pgen;
......
......@@ -141,9 +141,9 @@ private:
// get a map of mcID --> cluster
// input: cluster_collections --> list of handles to all cluster collections
std::map<eicd::Index, eicd::ConstCluster>
std::map<eicd::Index, eicd::Cluster>
indexedClusters(const std::vector<DataHandle<eicd::ClusterCollection>*>& cluster_collections) const {
std::map<eicd::Index, eicd::ConstCluster> matched = {};
std::map<eicd::Index, eicd::Cluster> matched = {};
for (const auto& cluster_handle : cluster_collections) {
const auto& clusters = *(cluster_handle->get());
for (const auto& cluster : clusters) {
......@@ -174,7 +174,7 @@ private:
// reconstruct a neutral cluster
// (for now assuming the vertex is at (0,0,0))
eicd::ReconstructedParticle reconstruct_neutral(const eicd::ConstCluster& clus, const double mass,
eicd::ReconstructedParticle reconstruct_neutral(const eicd::Cluster& clus, const double mass,
const int32_t pid) const {
const float energy = clus.energy();
const float momentum = energy < mass ? 0 : std::sqrt(energy * energy - mass * mass);
......
......@@ -64,8 +64,8 @@ public:
Rndm::Numbers m_gaussDist;
private:
using RecPart = eicd::ReconstructedParticle;
using Assoc = eicd::MCRecoParticleAssociation;
using RecPart = eicd::MutableReconstructedParticle;
using Assoc = eicd::MutableMCRecoParticleAssociation;
using RecData = std::pair<RecPart, Assoc>;
public:
......@@ -346,7 +346,7 @@ private:
// all momentum smearing in EIC-smear for the far-forward region uses
// the same 2 relations for P and Pt smearing (B0, RP, OMD)
RecData smearMomentum(const edm4hep::ConstMCParticle& part) {
RecData smearMomentum(const edm4hep::MCParticle& part) {
const auto mom_ion = rotateLabToIonDirection(part.getMomentum());
const double p = std::hypot(mom_ion.x, mom_ion.y, mom_ion.z);
const double dp = (0.025 * p) * m_gaussDist();
......@@ -367,7 +367,7 @@ private:
// And build our 3-vector
const edm4hep::Vector3f psmear_ion{static_cast<float>(pxs), static_cast<float>(pys), static_cast<float>(pzs)};
const auto psmear = rotateIonToLabDirection(psmear_ion);
eicd::ReconstructedParticle rec_part;
eicd::MutableReconstructedParticle rec_part;
rec_part.setType(-1);
rec_part.setEnergy(std::hypot(ps, part.getMass()));
rec_part.setMomentum({psmear.x, psmear.y, psmear.z});
......
......@@ -155,8 +155,8 @@ public:
// get a map of mcID --> cluster
// input: cluster_collections --> list of handles to all cluster collections
std::map<eicd::Index, eicd::ConstCluster> indexedClusters(const eicd::ClusterCollection& clusters) const {
std::map<eicd::Index, eicd::ConstCluster> matched = {};
std::map<eicd::Index, eicd::Cluster> indexedClusters(const eicd::ClusterCollection& clusters) const {
std::map<eicd::Index, eicd::Cluster> matched = {};
for (const auto& cluster : clusters) {
if (msgLevel(MSG::VERBOSE)) {
verbose() << " --> Found cluster: " << cluster.id() << " with mcID " << cluster.mcID() << " and energy "
......
......@@ -89,7 +89,7 @@ public:
auto& mhits = *m_outputHitCollection.createAndPut();
// container
std::unordered_map<std::pair<int64_t, int64_t>, std::vector<eicd::ConstCalorimeterHit>, pair_hash> merged_hits;
std::unordered_map<std::pair<int64_t, int64_t>, std::vector<eicd::CalorimeterHit>, pair_hash> merged_hits;
for (const auto h : *m_inputHitCollection.get()) {
auto bins =
......@@ -100,7 +100,7 @@ public:
for (const auto& [bins, hits] : merged_hits) {
const auto ref = hits.front();
eicd::CalorimeterHit hit;
eicd::MutableCalorimeterHit hit;
hit.setCellID(ref.getCellID());
// TODO, we can do timing cut to reject noises
hit.setTime(ref.getTime());
......
......@@ -107,7 +107,7 @@ public:
auto& outputs = *m_outputHitCollection.createAndPut();
// find the hits that belong to the same group (for merging)
std::unordered_map<long long, std::vector<eicd::ConstCalorimeterHit>> merge_map;
std::unordered_map<long long, std::vector<eicd::CalorimeterHit>> merge_map;
for (const auto& h : inputs) {
int64_t id = h.getCellID() & id_mask;
// use the reference field position
......
......@@ -41,7 +41,7 @@ using namespace Gaudi::Units;
namespace {
using CaloHit = eicd::ConstCalorimeterHit;
using CaloHit = eicd::CalorimeterHit;
using CaloHitCollection = eicd::CalorimeterHitCollection;
// helper functions to get distance between hits
......@@ -327,7 +327,7 @@ private:
}
return;
} else if (maxima.size() == 1) {
eicd::ProtoCluster pcl;
eicd::MutableProtoCluster pcl;
for (auto& [idx, hit] : group) {
pcl.addToHits(hit);
pcl.addToWeights(1.);
......@@ -342,7 +342,7 @@ private:
// split between maxima
// TODO, here we can implement iterations with profile, or even ML for better splits
std::vector<double> weights(maxima.size(), 1.);
std::vector<eicd::ProtoCluster> pcls;
std::vector<eicd::MutableProtoCluster> pcls;
for (size_t k = 0; k < maxima.size(); ++k) {
pcls.push_back({});
}
......
......@@ -139,8 +139,8 @@ public:
}
private:
eicd::Cluster reconstruct(const eicd::ConstProtoCluster& pcl) const {
eicd::Cluster cl;
eicd::MutableCluster reconstruct(const eicd::ProtoCluster& pcl) const {
eicd::MutableCluster cl;
cl.setNhits(pcl.hits_size());
// no hits
......
......@@ -252,7 +252,7 @@ public:
//----- end RP reconstruction code ------
eicd::ReconstructedParticle rpTrack;
eicd::MutableReconstructedParticle rpTrack;
rpTrack.setType(0);
rpTrack.setMomentum({prec});
rpTrack.setEnergy(std::hypot(eicd::magnitude(rpTrack.getMomentum()), .938272));
......
......@@ -159,7 +159,7 @@ public:
//----- end RP reconstruction code ------
eicd::ReconstructedParticle rpTrack;
eicd::MutableReconstructedParticle rpTrack;
rpTrack.setType(0);
rpTrack.setMomentum({prec});
rpTrack.setEnergy(std::hypot(eicd::magnitude(rpTrack.getMomentum()), .938272));
......
......@@ -131,11 +131,11 @@ public:
private:
template <typename T> static inline T pow2(const T& x) { return x * x; }
std::vector<eicd::Cluster> reconstruct_cluster_layers(const eicd::ConstProtoCluster& pcl) const {
std::vector<eicd::Cluster> reconstruct_cluster_layers(const eicd::ProtoCluster& pcl) const {
const auto& hits = pcl.getHits();
const auto& weights = pcl.getWeights();
// using map to have hits sorted by layer
std::map<int, std::vector<std::pair<eicd::ConstCalorimeterHit, float>>> layer_map;
std::map<int, std::vector<std::pair<eicd::CalorimeterHit, float>>> layer_map;
for (unsigned i = 0; i < hits.size(); ++i) {
const auto hit = hits[i];
auto lid = hit.getLayer();
......@@ -154,8 +154,8 @@ private:
return cl_layers;
}
eicd::Cluster reconstruct_layer(const std::vector<std::pair<eicd::ConstCalorimeterHit, float>>& hits) const {
eicd::Cluster layer;
eicd::Cluster reconstruct_layer(const std::vector<std::pair<eicd::CalorimeterHit, float>>& hits) const {
eicd::MutableCluster layer;
layer.setType(ClusterType::kClusterSlice);
// Calculate averages
double energy;
......@@ -193,8 +193,8 @@ private:
return layer;
}
eicd::Cluster reconstruct_cluster(const eicd::ConstProtoCluster& pcl) {
eicd::Cluster cluster;
eicd::MutableCluster reconstruct_cluster(const eicd::ProtoCluster& pcl) {
eicd::MutableCluster cluster;
const auto& hits = pcl.getHits();
const auto& weights = pcl.getWeights();
......
......@@ -74,7 +74,7 @@ namespace Jug::Reco {
auto& mhits = *m_outputHitCollection.createAndPut();
// group the hits by layer
std::vector<std::vector<eicd::ConstCalorimeterHit>> layer_hits;
std::vector<std::vector<eicd::CalorimeterHit>> layer_hits;
layer_hits.resize(m_nLayers);
for (const auto& h : hits) {
auto k = h.getLayer();
......@@ -86,7 +86,7 @@ namespace Jug::Reco {
// sort by energy
for (auto &layer : layer_hits) {
std::sort(layer.begin(), layer.end(),
[] (const eicd::ConstCalorimeterHit &h1, const eicd::ConstCalorimeterHit &h2) {
[] (const eicd::CalorimeterHit &h1, const eicd::CalorimeterHit &h2) {
return h1.getEnergy() > h2.getEnergy();
});
}
......
......@@ -131,7 +131,7 @@ public:
// group neighboring hits
std::vector<bool> visits(hits.size(), false);
std::vector<std::vector<std::pair<uint32_t, eicd::ConstCalorimeterHit>>> groups;
std::vector<std::vector<std::pair<uint32_t, eicd::CalorimeterHit>>> groups;
for (size_t i = 0; i < hits.size(); ++i) {
if (msgLevel(MSG::DEBUG)) {
debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})", i + 1,
......@@ -180,7 +180,7 @@ private:
template <typename T> static inline T pow2(const T& x) { return x * x; }
// helper function to group hits
bool is_neighbor(const eicd::ConstCalorimeterHit& h1, const eicd::ConstCalorimeterHit& h2) const {
bool is_neighbor(const eicd::CalorimeterHit& h1, const eicd::CalorimeterHit& h2) const {
// different sectors, simple distance check
if (h1.getSector() != h2.getSector()) {
return std::sqrt(pow2(h1.getPosition().x - h2.getPosition().x) + pow2(h1.getPosition().y - h2.getPosition().y) +
......@@ -204,7 +204,7 @@ private:
}
// grouping function with Depth-First Search
void dfs_group(std::vector<std::pair<uint32_t, eicd::ConstCalorimeterHit>>& group, int idx,
void dfs_group(std::vector<std::pair<uint32_t, eicd::CalorimeterHit>>& group, int idx,
const eicd::CalorimeterHitCollection& hits, std::vector<bool>& visits) const {
// not a qualified hit to participate in clustering, stop here
if (hits[idx].getEnergy() < minClusterHitEdep) {
......
......@@ -93,13 +93,13 @@ namespace Jug::Reco {
// mcHits = m_inputMC->get();
//}
std::vector<std::pair<uint32_t, eicd::ConstCalorimeterHit>> the_hits;
std::vector<std::pair<uint32_t, eicd::ConstCalorimeterHit>> remaining_hits;
std::vector<std::pair<uint32_t, eicd::CalorimeterHit>> the_hits;
std::vector<std::pair<uint32_t, eicd::CalorimeterHit>> remaining_hits;
double max_dist = m_maxDistance.value() / mm;
double min_energy = m_minModuleEdep.value() / GeV;
eicd::ConstCalorimeterHit ref_hit;
eicd::CalorimeterHit ref_hit;
bool have_ref = false;
// Collect all our hits, and get the highest energy hit
{
......@@ -120,7 +120,7 @@ namespace Jug::Reco {
while (have_ref && ref_hit.getEnergy() > min_energy) {
std::vector<std::pair<uint32_t, eicd::ConstCalorimeterHit>> cluster_hits;
std::vector<std::pair<uint32_t, eicd::CalorimeterHit>> cluster_hits;
for (const auto& [idx, h] : the_hits) {
if (eicd::magnitude(h.getPosition() - ref_hit.getPosition()) < max_dist) {
......@@ -132,7 +132,7 @@ namespace Jug::Reco {
double total_energy = std::accumulate(
std::begin(cluster_hits), std::end(cluster_hits), 0.0,
[](double t, const std::pair<uint32_t, eicd::ConstCalorimeterHit>& h1) { return (t + h1.second.getEnergy()); });
[](double t, const std::pair<uint32_t, eicd::CalorimeterHit>& h1) { return (t + h1.second.getEnergy()); });
if (msgLevel(MSG::DEBUG)) {
debug() << " total_energy = " << total_energy << endmsg;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment