diff --git a/src/algorithms/digi/SiliconTrackerDigi.cc b/src/algorithms/digi/SiliconTrackerDigi.cc index 5833f5ae17..b5f21216eb 100644 --- a/src/algorithms/digi/SiliconTrackerDigi.cc +++ b/src/algorithms/digi/SiliconTrackerDigi.cc @@ -20,10 +20,7 @@ namespace eicrecon { -void SiliconTrackerDigi::init(std::shared_ptr& logger) { - // set logger - m_log = logger; - +void SiliconTrackerDigi::init() { // Create random gauss function m_gauss = [&](){ return m_random.Gaus(0, m_cfg.timeResolution); @@ -50,20 +47,20 @@ void SiliconTrackerDigi::process( double result_time = sim_hit.getTime() + time_smearing; auto hit_time_stamp = (std::int32_t) (result_time * 1e3); - m_log->debug("--------------------"); - m_log->debug("Hit cellID = {}", sim_hit.getCellID()); - m_log->debug(" position = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z); - m_log->debug(" xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y)); - m_log->debug(" momentum = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z); - m_log->debug(" edep = {:.2f}", sim_hit.getEDep()); - m_log->debug(" time = {:.4f}[ns]", sim_hit.getTime()); - m_log->debug(" particle time = {}[ns]", sim_hit.getMCParticle().getTime()); - m_log->debug(" time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time); - m_log->debug(" hit_time_stamp: {} [~ps]", hit_time_stamp); + debug("--------------------"); + debug("Hit cellID = {}", sim_hit.getCellID()); + debug(" position = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getPosition().x, sim_hit.getPosition().y, sim_hit.getPosition().z); + debug(" xy_radius = {:.2f}", std::hypot(sim_hit.getPosition().x, sim_hit.getPosition().y)); + debug(" momentum = ({:.2f}, {:.2f}, {:.2f})", sim_hit.getMomentum().x, sim_hit.getMomentum().y, sim_hit.getMomentum().z); + debug(" edep = {:.2f}", sim_hit.getEDep()); + debug(" time = {:.4f}[ns]", sim_hit.getTime()); + debug(" particle time = {}[ns]", sim_hit.getMCParticle().getTime()); + debug(" time smearing: {:.4f}, resulting time = {:.4f} [ns]", time_smearing, result_time); + debug(" hit_time_stamp: {} [~ps]", hit_time_stamp); if (sim_hit.getEDep() < m_cfg.threshold) { - m_log->debug(" edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV); + debug(" edep is below threshold of {:.2f} [keV]", m_cfg.threshold / dd4hep::keV); continue; } @@ -77,7 +74,7 @@ void SiliconTrackerDigi::process( } else { // There is previous values in the cell auto& hit = cell_hit_map[sim_hit.getCellID()]; - m_log->debug(" Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp()); + debug(" Hit already exists in cell ID={}, prev. hit time: {}", sim_hit.getCellID(), hit.getTimeStamp()); // keep earliest time for hit auto time_stamp = hit.getTimeStamp(); diff --git a/src/algorithms/digi/SiliconTrackerDigi.h b/src/algorithms/digi/SiliconTrackerDigi.h index 1280783a9f..52afcd44e1 100644 --- a/src/algorithms/digi/SiliconTrackerDigi.h +++ b/src/algorithms/digi/SiliconTrackerDigi.h @@ -5,12 +5,10 @@ #include #include -#include #include +#include #include -#include #include -#include #include #include @@ -19,45 +17,35 @@ namespace eicrecon { - using SiliconTrackerDigiAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::SimTrackerHitCollection - >, - algorithms::Output< - edm4eic::RawTrackerHitCollection, - edm4eic::MCRecoTrackerHitAssociationCollection - > - >; - - class SiliconTrackerDigi - : public SiliconTrackerDigiAlgorithm, - public WithPodConfig { - - public: - SiliconTrackerDigi(std::string_view name) - : SiliconTrackerDigiAlgorithm{name, - {"inputHitCollection"}, - {"outputRawHitCollection","outputHitAssociations"}, - "Apply threshold, digitize within ADC range, " - "convert time with smearing resolution."} {} +using SiliconTrackerDigiAlgorithm = + algorithms::Algorithm, + algorithms::Output>; - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; +class SiliconTrackerDigi : public SiliconTrackerDigiAlgorithm, + public WithPodConfig { - private: - /** algorithm logger */ - std::shared_ptr m_log; +public: + SiliconTrackerDigi(std::string_view name) + : SiliconTrackerDigiAlgorithm{name, + {"inputHitCollection"}, + {"outputRawHitCollection", "outputHitAssociations"}, + "Apply threshold, digitize within ADC range, " + "convert time with smearing resolution."} {} - /** Random number generation*/ - TRandomMixMax m_random; - std::function m_gauss; + void init() final; + void process(const Input&, const Output&) const final; - // FIXME replace with standard random engine - //std::default_random_engine generator; // TODO: need something more appropriate here - //std::normal_distribution m_normDist; // defaults to mean=0, sigma=1 +private: + /** Random number generation*/ + TRandomMixMax m_random; + std::function m_gauss; - //algorithms::Generator m_rng = algorithms::RandomSvc::instance().generator(); + // FIXME replace with standard random engine + // std::default_random_engine generator; // TODO: need something more appropriate here + // std::normal_distribution m_normDist; // defaults to mean=0, sigma=1 - }; + // algorithms::Generator m_rng = algorithms::RandomSvc::instance().generator(); +}; -} // eicrecon +} // namespace eicrecon diff --git a/src/algorithms/fardetectors/FarDetectorLinearTracking.cc b/src/algorithms/fardetectors/FarDetectorLinearTracking.cc index e323ff6db9..5d93d42bce 100644 --- a/src/algorithms/fardetectors/FarDetectorLinearTracking.cc +++ b/src/algorithms/fardetectors/FarDetectorLinearTracking.cc @@ -25,9 +25,7 @@ namespace eicrecon { - void FarDetectorLinearTracking::init(std::shared_ptr& logger) { - - m_log = logger; + void FarDetectorLinearTracking::init() { // For changing how strongly each layer hit is in contributing to the fit m_layerWeights = Eigen::VectorXd::Constant(m_cfg.n_layer,1); @@ -49,7 +47,7 @@ namespace eicrecon { // Check the number of input collections is correct int nCollections = inputhits.size(); if(nCollections!=m_cfg.n_layer){ - m_log->error("Wrong number of input collections passed to algorithm"); + error("Wrong number of input collections passed to algorithm"); return; } @@ -58,7 +56,7 @@ namespace eicrecon { // TODO - Implement more sensible solution for(const auto& layerHits: inputhits){ if((*layerHits).size()>m_cfg.layer_hits_max){ - m_log->info("Too many hits in layer"); + info("Too many hits in layer"); return; } } @@ -160,9 +158,9 @@ namespace eicrecon { double angle = std::acos(hitDiff.dot(m_optimumDirection)); - m_log->debug("Vector: x={}, y={}, z={}",hitDiff.x(),hitDiff.y(),hitDiff.z()); - m_log->debug("Optimum: x={}, y={}, z={}",m_optimumDirection.x(),m_optimumDirection.y(),m_optimumDirection.z()); - m_log->debug("Angle: {}, Tolerance {}",angle,m_cfg.step_angle_tolerance); + debug("Vector: x={}, y={}, z={}",hitDiff.x(),hitDiff.y(),hitDiff.z()); + debug("Optimum: x={}, y={}, z={}",m_optimumDirection.x(),m_optimumDirection.y(),m_optimumDirection.z()); + debug("Angle: {}, Tolerance {}",angle,m_cfg.step_angle_tolerance); if(angle>m_cfg.step_angle_tolerance) return false; diff --git a/src/algorithms/fardetectors/FarDetectorLinearTracking.h b/src/algorithms/fardetectors/FarDetectorLinearTracking.h index 5d4013373c..d9caea23f5 100644 --- a/src/algorithms/fardetectors/FarDetectorLinearTracking.h +++ b/src/algorithms/fardetectors/FarDetectorLinearTracking.h @@ -3,65 +3,54 @@ #pragma once +#include #include #include #include #include #include -#include #include #include #include -#include -#include + #include "FarDetectorLinearTrackingConfig.h" namespace eicrecon { - using FarDetectorLinearTrackingAlgorithm = algorithms::Algorithm< - algorithms::Input< - std::vector - >, - algorithms::Output< - edm4eic::TrackSegmentCollection - > - >; - - class FarDetectorLinearTracking - : public FarDetectorLinearTrackingAlgorithm, - public WithPodConfig { - - public: - FarDetectorLinearTracking(std::string_view name) - : FarDetectorLinearTrackingAlgorithm{name, - {"inputHitCollections"}, - {"outputTrackSegments"}, - "Fit track segments from hits in the tracker layers"} {} +using FarDetectorLinearTrackingAlgorithm = + algorithms::Algorithm>, + algorithms::Output>; - /** One time initialization **/ - void init(std::shared_ptr& logger); +class FarDetectorLinearTracking : public FarDetectorLinearTrackingAlgorithm, + public WithPodConfig { - /** Event by event processing **/ - void process(const Input&, const Output&) const final; +public: + FarDetectorLinearTracking(std::string_view name) + : FarDetectorLinearTrackingAlgorithm{name, + {"inputHitCollections"}, + {"outputTrackSegments"}, + "Fit track segments from hits in the tracker layers"} {} - private: - std::shared_ptr m_log; + /** One time initialization **/ + void init() final; - Eigen::VectorXd m_layerWeights; + /** Event by event processing **/ + void process(const Input&, const Output&) const final; - Eigen::Vector3d m_optimumDirection; +private: + Eigen::VectorXd m_layerWeights; - void buildMatrixRecursive(int level, - Eigen::MatrixXd* hitMatrix, - const std::vector>& hits, - gsl::not_null outputTracks) const; + Eigen::Vector3d m_optimumDirection; - void checkHitCombination(Eigen::MatrixXd* hitMatrix, - gsl::not_null outputTracks) const; + void + buildMatrixRecursive(int level, Eigen::MatrixXd* hitMatrix, + const std::vector>& hits, + gsl::not_null outputTracks) const; - bool checkHitPair(const Eigen::Vector3d& hit1, - const Eigen::Vector3d& hit2) const; + void checkHitCombination(Eigen::MatrixXd* hitMatrix, + gsl::not_null outputTracks) const; - }; + bool checkHitPair(const Eigen::Vector3d& hit1, const Eigen::Vector3d& hit2) const; +}; -} // eicrecon +} // namespace eicrecon diff --git a/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc b/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc index 287dea71e4..b12a61b738 100644 --- a/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc +++ b/src/algorithms/fardetectors/FarDetectorTrackerCluster.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -15,192 +16,186 @@ #include #include #include -#include #include +#include #include "algorithms/fardetectors/FarDetectorTrackerCluster.h" #include "algorithms/fardetectors/FarDetectorTrackerClusterConfig.h" namespace eicrecon { - void FarDetectorTrackerCluster::init(std::shared_ptr& log) { +void FarDetectorTrackerCluster::init() { - m_log = log; - m_detector = algorithms::GeoSvc::instance().detector(); - m_cellid_converter = algorithms::GeoSvc::instance().cellIDPositionConverter(); + m_detector = algorithms::GeoSvc::instance().detector(); + m_cellid_converter = algorithms::GeoSvc::instance().cellIDPositionConverter(); - if (m_cfg.readout.empty()) { - throw JException("Readout is empty"); + if (m_cfg.readout.empty()) { + throw JException("Readout is empty"); + } + try { + m_seg = m_detector->readout(m_cfg.readout).segmentation(); + m_id_dec = m_detector->readout(m_cfg.readout).idSpec().decoder(); + if (!m_cfg.x_field.empty()) { + m_x_idx = m_id_dec->index(m_cfg.x_field); + debug("Find layer field {}, index = {}", m_cfg.x_field, m_x_idx); } - try { - m_seg = m_detector->readout(m_cfg.readout).segmentation(); - m_id_dec = m_detector->readout(m_cfg.readout).idSpec().decoder(); - if (!m_cfg.x_field.empty()) { - m_x_idx = m_id_dec->index(m_cfg.x_field); - m_log->debug("Find layer field {}, index = {}", m_cfg.x_field, m_x_idx); - } - if (!m_cfg.y_field.empty()) { - m_y_idx = m_id_dec->index(m_cfg.y_field); - m_log->debug("Find layer field {}, index = {}", m_cfg.y_field, m_y_idx); - } - } catch (...) { - m_log->error("Failed to load ID decoder for {}", m_cfg.readout); - throw JException("Failed to load ID decoder"); + if (!m_cfg.y_field.empty()) { + m_y_idx = m_id_dec->index(m_cfg.y_field); + debug("Find layer field {}, index = {}", m_cfg.y_field, m_y_idx); } - + } catch (...) { + error("Failed to load ID decoder for {}", m_cfg.readout); + throw JException("Failed to load ID decoder"); } +} - void FarDetectorTrackerCluster::process( - const FarDetectorTrackerCluster::Input& input, - const FarDetectorTrackerCluster::Output& output) const { - - const auto [inputHitsCollections] = input; - auto [outputClustersCollection] = output; +void FarDetectorTrackerCluster::process(const FarDetectorTrackerCluster::Input& input, + const FarDetectorTrackerCluster::Output& output) const { - //Loop over input and output collections - Any collection should only contain hits from a single surface - for(size_t i=0; isize() == 0) continue; - auto outputClusters = outputClustersCollection[i]; + const auto [inputHitsCollections] = input; + auto [outputClustersCollection] = output; - // Make clusters - auto clusters = ClusterHits(*inputHits); + // Loop over input and output collections - Any collection should only contain hits from a single + // surface + for (size_t i = 0; i < inputHitsCollections.size(); i++) { + auto inputHits = inputHitsCollections[i]; + if (inputHits->size() == 0) + continue; + auto outputClusters = outputClustersCollection[i]; - // Create TrackerHits from 2D cluster positions - ConvertClusters(clusters, *outputClusters); + // Make clusters + auto clusters = ClusterHits(*inputHits); - } + // Create TrackerHits from 2D cluster positions + ConvertClusters(clusters, *outputClusters); } +} - // Create vector of FDTrackerCluster from list of hits - std::vector FarDetectorTrackerCluster::ClusterHits(const edm4eic::RawTrackerHitCollection& inputHits) const { - - std::vector clusters; - - ROOT::VecOps::RVec id; - ROOT::VecOps::RVec x; - ROOT::VecOps::RVec y; - ROOT::VecOps::RVec e; - ROOT::VecOps::RVec t; - - // Gather detector id positions - for(const auto& hit: inputHits){ - auto cellID = hit.getCellID(); - id.push_back(cellID); - x.push_back (m_id_dec->get( cellID, m_x_idx )); - y.push_back (m_id_dec->get( cellID, m_y_idx )); - e.push_back (hit.getCharge()); - t.push_back (hit.getTimeStamp()); - } - - // Set up clustering variables - ROOT::VecOps::RVec available(id.size(), 1); - auto indices = Enumerate(id); - - // Loop while there are unclustered hits - while(ROOT::VecOps::Any(available)){ - - dd4hep::Position localPos = {0,0,0}; - float weightSum = 0; - - float esum = 0; - float t0 = 0; - float tError = 0; - auto maxIndex = ROOT::VecOps::ArgMax(e*available); - - available[maxIndex] = 0; - - ROOT::VecOps::RVec clusterList = {maxIndex}; - ROOT::VecOps::RVec clusterT; - std::vector clusterHits; +// Create vector of FDTrackerCluster from list of hits +std::vector +FarDetectorTrackerCluster::ClusterHits(const edm4eic::RawTrackerHitCollection& inputHits) const { + + std::vector clusters; + + ROOT::VecOps::RVec id; + ROOT::VecOps::RVec x; + ROOT::VecOps::RVec y; + ROOT::VecOps::RVec e; + ROOT::VecOps::RVec t; + + // Gather detector id positions + for (const auto& hit : inputHits) { + auto cellID = hit.getCellID(); + id.push_back(cellID); + x.push_back(m_id_dec->get(cellID, m_x_idx)); + y.push_back(m_id_dec->get(cellID, m_y_idx)); + e.push_back(hit.getCharge()); + t.push_back(hit.getTimeStamp()); + } - // Loop over hits, adding neighbouring hits as relevant - while(clusterList.size()){ + // Set up clustering variables + ROOT::VecOps::RVec available(id.size(), 1); + auto indices = Enumerate(id); - // Takes first remaining hit in cluster list - auto index = clusterList[0]; + // Loop while there are unclustered hits + while (ROOT::VecOps::Any(available)) { - // Finds neighbours of cluster within time limit - auto filter = available*(abs(x-x[index])<=1)*(abs(y-y[index])<=1)*(abs(t-t[index]) clusterList = {maxIndex}; + ROOT::VecOps::RVec clusterT; + std::vector clusterHits; - // Adds raw hit to TrackerHit contribution - clusterHits.push_back((inputHits)[index].getObjectID()); + // Loop over hits, adding neighbouring hits as relevant + while (clusterList.size()) { - // Energy - auto hitE = e[index]; - esum += hitE; - // TODO - See if now a single detector element is expected a better function is available. - auto pos = m_seg->position(id[index]); + // Takes first remaining hit in cluster list + auto index = clusterList[0]; - //Weighted position - float weight = hitE; // TODO - Calculate appropriate weighting based on sensor charge sharing - weightSum += weight; - localPos += pos*weight; + // Finds neighbours of cluster within time limit + auto filter = available * (abs(x - x[index]) <= 1) * (abs(y - y[index]) <= 1) * + (abs(t - t[index]) < m_cfg.hit_time_limit); - //Time - clusterT.push_back(t[index]); + // Adds the found hits to the cluster + clusterList = Concatenate(clusterList, indices[filter]); - } + // Removes the found hits from the list of still available hits + available = available * (!filter); - // Finalise position - localPos/=weightSum; + // Removes current hit from remaining found cluster hits + clusterList.erase(clusterList.begin()); - // Finalise time - t0 = Mean(clusterT); - tError = StdDev(clusterT); // TODO fold detector timing resolution into error + // Adds raw hit to TrackerHit contribution + clusterHits.push_back((inputHits)[index].getObjectID()); - // Create cluster - clusters.push_back(FDTrackerCluster{ - .cellID = id[maxIndex], - .x = localPos.x(), - .y = localPos.y(), - .energy = esum, - .time = t0, - .timeError = tError, - .rawHits = clusterHits - }); + // Energy + auto hitE = e[index]; + esum += hitE; + // TODO - See if now a single detector element is expected a better function is available. + auto pos = m_seg->position(id[index]); + // Weighted position + float weight = hitE; // TODO - Calculate appropriate weighting based on sensor charge sharing + weightSum += weight; + localPos += pos * weight; + // Time + clusterT.push_back(t[index]); } - return clusters; - + // Finalise position + localPos /= weightSum; + + // Finalise time + t0 = Mean(clusterT); + tError = StdDev(clusterT); // TODO fold detector timing resolution into error + + // Create cluster + clusters.push_back(FDTrackerCluster{.cellID = id[maxIndex], + .x = localPos.x(), + .y = localPos.y(), + .energy = esum, + .time = t0, + .timeError = tError, + .rawHits = clusterHits}); } - // Convert to global coordinates and create TrackerHits - void FarDetectorTrackerCluster::ConvertClusters(const std::vector& clusters, edm4hep::TrackerHitCollection& outputClusters) const { + return clusters; +} - // Get context of first hit - const dd4hep::VolumeManagerContext* context = m_cellid_converter->findContext(clusters[0].cellID); +// Convert to global coordinates and create TrackerHits +void FarDetectorTrackerCluster::ConvertClusters( + const std::vector& clusters, + edm4hep::TrackerHitCollection& outputClusters) const { - for(auto cluster: clusters){ - auto hitPos = outputClusters.create(); + // Get context of first hit + const dd4hep::VolumeManagerContext* context = m_cellid_converter->findContext(clusters[0].cellID); - auto globalPos = context->localToWorld({cluster.x,cluster.y,0}); + for (auto cluster : clusters) { + auto hitPos = outputClusters.create(); - // Set cluster members - hitPos.setCellID (cluster.cellID); - hitPos.setPosition(edm4hep::Vector3d(globalPos.x()/dd4hep::mm,globalPos.y()/dd4hep::mm,globalPos.z()/dd4hep::mm)); - hitPos.setEDep (cluster.energy); - hitPos.setTime (cluster.time); + auto globalPos = context->localToWorld({cluster.x, cluster.y, 0}); - // Add raw hits to cluster - for(auto hit: cluster.rawHits){ - hitPos.addToRawHits(hit); - } + // Set cluster members + hitPos.setCellID(cluster.cellID); + hitPos.setPosition(edm4hep::Vector3d(globalPos.x() / dd4hep::mm, globalPos.y() / dd4hep::mm, + globalPos.z() / dd4hep::mm)); + hitPos.setEDep(cluster.energy); + hitPos.setTime(cluster.time); + // Add raw hits to cluster + for (auto hit : cluster.rawHits) { + hitPos.addToRawHits(hit); } - } - - } + +} // namespace eicrecon diff --git a/src/algorithms/fardetectors/FarDetectorTrackerCluster.h b/src/algorithms/fardetectors/FarDetectorTrackerCluster.h index 381eb97886..6ac60e5067 100644 --- a/src/algorithms/fardetectors/FarDetectorTrackerCluster.h +++ b/src/algorithms/fardetectors/FarDetectorTrackerCluster.h @@ -11,8 +11,6 @@ #include #include #include -#include -#include #include #include #include @@ -32,48 +30,41 @@ struct FDTrackerCluster { }; namespace eicrecon { - using FarDetectorTrackerClusterAlgorithm = algorithms::Algorithm< - algorithms::Input< - std::vector - >, - algorithms::Output< - std::vector - > - >; +using FarDetectorTrackerClusterAlgorithm = + algorithms::Algorithm>, + algorithms::Output>>; - class FarDetectorTrackerCluster - : public FarDetectorTrackerClusterAlgorithm, - public WithPodConfig { +class FarDetectorTrackerCluster : public FarDetectorTrackerClusterAlgorithm, + public WithPodConfig { - public: - FarDetectorTrackerCluster(std::string_view name) +public: + FarDetectorTrackerCluster(std::string_view name) : FarDetectorTrackerClusterAlgorithm{name, - {"inputHitCollection"}, - {"outputClusterPositionCollection"}, - "Simple weighted clustering of hits by x-y component of single detector element segmentation"} {} + {"inputHitCollection"}, + {"outputClusterPositionCollection"}, + "Simple weighted clustering of hits by x-y component of " + "single detector element segmentation"} {} - /** One time initialization **/ - void init(std::shared_ptr& logger); + /** One time initialization **/ + void init() final; - /** Event by event processing **/ - void process(const Input&, const Output&) const final; + /** Event by event processing **/ + void process(const Input&, const Output&) const final; - /** Cluster hits **/ - std::vector ClusterHits(const edm4eic::RawTrackerHitCollection&) const; + /** Cluster hits **/ + std::vector ClusterHits(const edm4eic::RawTrackerHitCollection&) const; - /** Convert clusters to TrackerHits **/ - void ConvertClusters(const std::vector&, edm4hep::TrackerHitCollection&) const; + /** Convert clusters to TrackerHits **/ + void ConvertClusters(const std::vector&, edm4hep::TrackerHitCollection&) const; - private: - const dd4hep::Detector* m_detector{nullptr}; - const dd4hep::BitFieldCoder* m_id_dec{nullptr}; - std::shared_ptr m_log; - const dd4hep::rec::CellIDPositionConverter* m_cellid_converter{nullptr}; - dd4hep::Segmentation m_seg; +private: + const dd4hep::Detector* m_detector{nullptr}; + const dd4hep::BitFieldCoder* m_id_dec{nullptr}; + const dd4hep::rec::CellIDPositionConverter* m_cellid_converter{nullptr}; + dd4hep::Segmentation m_seg; - int m_x_idx{0}; - int m_y_idx{0}; - - }; + int m_x_idx{0}; + int m_y_idx{0}; +}; -} // eicrecon +} // namespace eicrecon diff --git a/src/algorithms/meta/CollectionCollector.h b/src/algorithms/meta/CollectionCollector.h index accf942f42..f6806366b8 100644 --- a/src/algorithms/meta/CollectionCollector.h +++ b/src/algorithms/meta/CollectionCollector.h @@ -28,9 +28,7 @@ namespace eicrecon { "Merge content of collections into one subset collection" }{} - void init(std::shared_ptr& logger){ // set logger - m_log = logger; - }; + void init() final { }; void process(const typename CollectionCollector::Input& input, const typename CollectionCollector::Output& output) const final{ @@ -46,8 +44,5 @@ namespace eicrecon { } } - private: - std::shared_ptr m_log; // logger - }; } // eicrecon diff --git a/src/algorithms/onnx/InclusiveKinematicsML.cc b/src/algorithms/onnx/InclusiveKinematicsML.cc index 6ef2e72a40..935efadb77 100644 --- a/src/algorithms/onnx/InclusiveKinematicsML.cc +++ b/src/algorithms/onnx/InclusiveKinematicsML.cc @@ -30,9 +30,7 @@ namespace eicrecon { return tensor; } - void InclusiveKinematicsML::init(std::shared_ptr& logger) { - m_log = logger; - + void InclusiveKinematicsML::init() { // onnxruntime setup Ort::Env env(ORT_LOGGING_LEVEL_WARNING, "inclusive-kinematics-ml"); Ort::SessionOptions session_options; @@ -41,19 +39,19 @@ namespace eicrecon { // print name/shape of inputs Ort::AllocatorWithDefaultOptions allocator; - m_log->debug("Input Node Name/Shape:"); + debug("Input Node Name/Shape:"); for (std::size_t i = 0; i < m_session.GetInputCount(); i++) { m_input_names.emplace_back(m_session.GetInputNameAllocated(i, allocator).get()); m_input_shapes.emplace_back(m_session.GetInputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); - m_log->debug("\t{} : {}", m_input_names.at(i), print_shape(m_input_shapes.at(i))); + debug("\t{} : {}", m_input_names.at(i), print_shape(m_input_shapes.at(i))); } // print name/shape of outputs - m_log->debug("Output Node Name/Shape:"); + debug("Output Node Name/Shape:"); for (std::size_t i = 0; i < m_session.GetOutputCount(); i++) { m_output_names.emplace_back(m_session.GetOutputNameAllocated(i, allocator).get()); m_output_shapes.emplace_back(m_session.GetOutputTypeInfo(i).GetTensorTypeAndShapeInfo().GetShape()); - m_log->debug("\t{} : {}", m_output_names.at(i), print_shape(m_output_shapes.at(i))); + debug("\t{} : {}", m_output_names.at(i), print_shape(m_output_shapes.at(i))); } // convert names to char* @@ -65,7 +63,7 @@ namespace eicrecon { [&](const std::string& str) { return str.c_str(); }); } catch(std::exception& e) { - m_log->error(e.what()); + error(e.what()); } } @@ -78,13 +76,13 @@ namespace eicrecon { // Require valid inputs if (electron->size() == 0 || da->size() == 0) { - m_log->debug("skipping because input collections have no entries"); + debug("skipping because input collections have no entries"); return; } // Assume model has 1 input nodes and 1 output node. if (m_input_names.size() != 1 || m_output_names.size() != 1) { - m_log->debug("skipping because model has incorrect input and output size"); + debug("skipping because model has incorrect input and output size"); return; } @@ -98,7 +96,7 @@ namespace eicrecon { // Double-check the dimensions of the input tensor if (! input_tensors[0].IsTensor() || input_tensors[0].GetTensorTypeAndShapeInfo().GetShape() != m_input_shapes.front()) { - m_log->debug("skipping because input tensor shape incorrect"); + debug("skipping because input tensor shape incorrect"); return; } @@ -109,7 +107,7 @@ namespace eicrecon { // Double-check the dimensions of the output tensors if (!output_tensors[0].IsTensor() || output_tensors.size() != m_output_names.size()) { - m_log->debug("skipping because output tensor shape incorrect"); + debug("skipping because output tensor shape incorrect"); return; } @@ -120,7 +118,7 @@ namespace eicrecon { kin.setX(x); } catch (const Ort::Exception& exception) { - m_log->error("error running model inference: {}", exception.what()); + error("error running model inference: {}", exception.what()); } } diff --git a/src/algorithms/onnx/InclusiveKinematicsML.h b/src/algorithms/onnx/InclusiveKinematicsML.h index e2d559dbf0..485fdb5d8b 100644 --- a/src/algorithms/onnx/InclusiveKinematicsML.h +++ b/src/algorithms/onnx/InclusiveKinematicsML.h @@ -4,11 +4,9 @@ #pragma once #include +#include #include #include -#include -#include -#include #include #include #include @@ -18,43 +16,35 @@ namespace eicrecon { - using InclusiveKinematicsMLAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4eic::InclusiveKinematicsCollection, - edm4eic::InclusiveKinematicsCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsML - : public InclusiveKinematicsMLAlgorithm, - public WithPodConfig { - - public: - InclusiveKinematicsML(std::string_view name) - : InclusiveKinematicsMLAlgorithm{name, - {"inclusiveKinematicsElectron", "inclusiveKinematicsDA"}, - {"inclusiveKinematicsML"}, - "Determine inclusive kinematics using combined ML method."} {} +using InclusiveKinematicsMLAlgorithm = + algorithms::Algorithm, + algorithms::Output>; - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; +class InclusiveKinematicsML : public InclusiveKinematicsMLAlgorithm, + public WithPodConfig { - private: - std::shared_ptr m_log; +public: + InclusiveKinematicsML(std::string_view name) + : InclusiveKinematicsMLAlgorithm{name, + {"inclusiveKinematicsElectron", "inclusiveKinematicsDA"}, + {"inclusiveKinematicsML"}, + "Determine inclusive kinematics using combined ML method."} { + } - mutable Ort::Session m_session{nullptr}; + void init() final; + void process(const Input&, const Output&) const final; - std::vector m_input_names; - std::vector m_input_names_char; - std::vector> m_input_shapes; +private: + mutable Ort::Session m_session{nullptr}; - std::vector m_output_names; - std::vector m_output_names_char; - std::vector> m_output_shapes; + std::vector m_input_names; + std::vector m_input_names_char; + std::vector> m_input_shapes; - }; + std::vector m_output_names; + std::vector m_output_names_char; + std::vector> m_output_shapes; +}; } // namespace eicrecon diff --git a/src/algorithms/pid/MergeTracks.cc b/src/algorithms/pid/MergeTracks.cc index 41087a91fb..847994332f 100644 --- a/src/algorithms/pid/MergeTracks.cc +++ b/src/algorithms/pid/MergeTracks.cc @@ -3,71 +3,64 @@ #include "MergeTracks.h" +#include +#include #include #include #include -#include -#include -#include +#include #include +#include #include #include namespace eicrecon { -void MergeTracks::init(std::shared_ptr& logger) -{ - m_log = logger; -} - -void MergeTracks::process( - const MergeTracks::Input& input, - const MergeTracks::Output& output) const { +void MergeTracks::process(const MergeTracks::Input& input, + const MergeTracks::Output& output) const { const auto [in_track_collections] = input; - auto [out_tracks] = output; + auto [out_tracks] = output; // logging - m_log->trace("{:=^70}"," call MergeTracks::AlgorithmProcess "); + trace("{:=^70}", " call MergeTracks::AlgorithmProcess "); // check that all input collections have the same size std::unordered_map in_track_collection_size_distribution; - for(const auto& in_track_collection : in_track_collections) { + for (const auto& in_track_collection : in_track_collections) { ++in_track_collection_size_distribution[in_track_collection->size()]; } if (in_track_collection_size_distribution.size() != 1) { std::vector in_track_collection_sizes; std::transform(in_track_collections.begin(), in_track_collections.end(), - std::back_inserter(in_track_collection_sizes), - [](const auto& in_track_collection) { return in_track_collection->size(); }); - m_log->error("cannot merge input track collections with different sizes {}", fmt::join(in_track_collection_sizes, ", ")); + std::back_inserter(in_track_collection_sizes), + [](const auto& in_track_collection) { return in_track_collection->size(); }); + error("cannot merge input track collections with different sizes {}", + fmt::join(in_track_collection_sizes, ", ")); return; } // loop over track collection elements std::size_t n_tracks = in_track_collection_size_distribution.begin()->first; - for(std::size_t i_track=0; i_trackcreate(); std::vector out_track_points; // loop over collections for this track, and add each track's points to `out_track_points` - for(const auto& in_track_collection : in_track_collections) { + for (const auto& in_track_collection : in_track_collections) { const auto& in_track = in_track_collection->at(i_track); - for(const auto& point : in_track.getPoints()) + for (const auto& point : in_track.getPoints()) out_track_points.push_back(point); } // sort all `out_track_points` by time - std::sort( - out_track_points.begin(), - out_track_points.end(), - [] (edm4eic::TrackPoint& a, edm4eic::TrackPoint& b) { return a.time < b.time; } - ); + std::sort(out_track_points.begin(), out_track_points.end(), + [](edm4eic::TrackPoint& a, edm4eic::TrackPoint& b) { return a.time < b.time; }); // add these sorted points to `out_track` - for(const auto& point : out_track_points) + for (const auto& point : out_track_points) out_track.addToPoints(point); /* FIXME: merge other members, such as `length` and `lengthError`; diff --git a/src/algorithms/pid/MergeTracks.h b/src/algorithms/pid/MergeTracks.h index abb5520e13..853ec1eb1c 100644 --- a/src/algorithms/pid/MergeTracks.h +++ b/src/algorithms/pid/MergeTracks.h @@ -12,36 +12,26 @@ // data model #include #include -#include -#include +#include +#include #include namespace eicrecon { - using MergeTracksAlgorithm = algorithms::Algorithm< - algorithms::Input< - std::vector - >, - algorithms::Output< - edm4eic::TrackSegmentCollection - > - >; +using MergeTracksAlgorithm = + algorithms::Algorithm>, + algorithms::Output>; - class MergeTracks - : public MergeTracksAlgorithm { +class MergeTracks : public MergeTracksAlgorithm { - public: - MergeTracks(std::string_view name) +public: + MergeTracks(std::string_view name) : MergeTracksAlgorithm{name, - {"inputTrackSegments"}, - {"outputTrackSegments"}, - "Effectively 'zip' the input track segments."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - - }; -} + {"inputTrackSegments"}, + {"outputTrackSegments"}, + "Effectively 'zip' the input track segments."} {} + + void init() final{}; + void process(const Input&, const Output&) const final; +}; +} // namespace eicrecon diff --git a/src/algorithms/reco/HadronicFinalState.cc b/src/algorithms/reco/HadronicFinalState.cc index 9cc84121b9..016770a85b 100644 --- a/src/algorithms/reco/HadronicFinalState.cc +++ b/src/algorithms/reco/HadronicFinalState.cc @@ -25,11 +25,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void HadronicFinalState::init(std::shared_ptr& logger) { - m_log = logger; + void HadronicFinalState::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -44,7 +43,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -58,7 +57,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -72,7 +71,7 @@ namespace eicrecon { // Get first scattered electron const auto ef_coll = find_first_scattered_electron(mcparts); if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); + debug("No truth scattered electron found"); return; } // Associate first scattered electron with reconstructed electrons @@ -87,7 +86,7 @@ namespace eicrecon { } } if (!(ef_assoc != rcassoc->end())) { - m_log->debug("Truth scattered electron not in reconstructed particles"); + debug("Truth scattered electron not in reconstructed particles"); return; } const auto ef_rc{ef_assoc->getRec()}; @@ -133,11 +132,11 @@ namespace eicrecon { // Sigma zero or negative if (sigma <= 0) { - m_log->debug("Sigma zero or negative"); + debug("Sigma zero or negative"); return; } - m_log->debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma()); + debug("sigma_h, pT_h, gamma_h = {},{},{}", hfs.getSigma(), hfs.getPT(), hfs.getGamma()); } diff --git a/src/algorithms/reco/HadronicFinalState.h b/src/algorithms/reco/HadronicFinalState.h index e4de62ebc9..6b25f48d7f 100644 --- a/src/algorithms/reco/HadronicFinalState.h +++ b/src/algorithms/reco/HadronicFinalState.h @@ -8,41 +8,30 @@ #include #include #include -#include -#include #include #include - namespace eicrecon { - using HadronicFinalStateAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::MCRecoParticleAssociationCollection - >, - algorithms::Output< - edm4eic::HadronicFinalStateCollection - > - >; - - class HadronicFinalState - : public HadronicFinalStateAlgorithm { - - public: - HadronicFinalState(std::string_view name) +using HadronicFinalStateAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class HadronicFinalState : public HadronicFinalStateAlgorithm { + +public: + HadronicFinalState(std::string_view name) : HadronicFinalStateAlgorithm{name, - {"MCParticles", "inputParticles", "inputAssociations"}, - {"hadronicFinalState"}, - "Calculate summed quantities of the hadronic final state."} {} + {"MCParticles", "inputParticles", "inputAssociations"}, + {"hadronicFinalState"}, + "Calculate summed quantities of the hadronic final state."} {} - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; + void init() final; + void process(const Input&, const Output&) const final; - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsDA.cc b/src/algorithms/reco/InclusiveKinematicsDA.cc index dd58f36a45..d925f07f01 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.cc +++ b/src/algorithms/reco/InclusiveKinematicsDA.cc @@ -23,8 +23,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsDA::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsDA::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { // m_log->debug("Unable to locate Particle Service. " @@ -42,7 +41,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -56,7 +55,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -83,7 +82,7 @@ namespace eicrecon { // Sigma zero or negative if (sigma_h <= 0) { - m_log->debug("Sigma zero or negative"); + debug("Sigma zero or negative"); return; } @@ -96,7 +95,7 @@ namespace eicrecon { auto kin = kinematics->create(x_da, Q2_da, W_da, y_da, nu_da); kin.setScat(kf); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsDA.h b/src/algorithms/reco/InclusiveKinematicsDA.h index c1d3f72468..ab6b71f7de 100644 --- a/src/algorithms/reco/InclusiveKinematicsDA.h +++ b/src/algorithms/reco/InclusiveKinematicsDA.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsDAAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsDA - : public InclusiveKinematicsDAAlgorithm { - - public: - InclusiveKinematicsDA(std::string_view name) - : InclusiveKinematicsDAAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using double-angle method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsDAAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsDA : public InclusiveKinematicsDAAlgorithm { + +public: + InclusiveKinematicsDA(std::string_view name) + : InclusiveKinematicsDAAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using double-angle method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.cc b/src/algorithms/reco/InclusiveKinematicsElectron.cc index 0a4781db91..8ac014f2ff 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.cc +++ b/src/algorithms/reco/InclusiveKinematicsElectron.cc @@ -23,11 +23,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsElectron::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsElectron::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -90,7 +89,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -104,7 +103,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -124,7 +123,7 @@ namespace eicrecon { // If no scattered electron was found if (electrons.size() == 0) { - m_log->debug("No scattered electron found"); + debug("No scattered electron found"); return; } @@ -140,7 +139,7 @@ namespace eicrecon { auto kin = kinematics->create(x, Q2, W, y, nu); kin.setScat(escat->at(0)); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsElectron.h b/src/algorithms/reco/InclusiveKinematicsElectron.h index e95cc55138..8b4864d14a 100644 --- a/src/algorithms/reco/InclusiveKinematicsElectron.h +++ b/src/algorithms/reco/InclusiveKinematicsElectron.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsElectronAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsElectron - : public InclusiveKinematicsElectronAlgorithm { - - public: - InclusiveKinematicsElectron(std::string_view name) - : InclusiveKinematicsElectronAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using electron method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsElectronAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsElectron : public InclusiveKinematicsElectronAlgorithm { + +public: + InclusiveKinematicsElectron(std::string_view name) + : InclusiveKinematicsElectronAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using electron method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsJB.cc b/src/algorithms/reco/InclusiveKinematicsJB.cc index 324093b1aa..48874d79bd 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.cc +++ b/src/algorithms/reco/InclusiveKinematicsJB.cc @@ -20,11 +20,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsJB::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsJB::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -39,7 +38,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -53,7 +52,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -71,7 +70,7 @@ namespace eicrecon { // Sigma zero or negative if (sigma_h <= 0) { - m_log->debug("Sigma zero or negative"); + debug("Sigma zero or negative"); return; } @@ -84,7 +83,7 @@ namespace eicrecon { auto kin = kinematics->create(x_jb, Q2_jb, W_jb, y_jb, nu_jb); kin.setScat(escat->at(0)); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsJB.h b/src/algorithms/reco/InclusiveKinematicsJB.h index 24526c04b2..bacc5a69f7 100644 --- a/src/algorithms/reco/InclusiveKinematicsJB.h +++ b/src/algorithms/reco/InclusiveKinematicsJB.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsJBAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsJB - : public InclusiveKinematicsJBAlgorithm { - - public: - InclusiveKinematicsJB(std::string_view name) - : InclusiveKinematicsJBAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using Jacquet-Blondel method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsJBAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsJB : public InclusiveKinematicsJBAlgorithm { + +public: + InclusiveKinematicsJB(std::string_view name) + : InclusiveKinematicsJBAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using Jacquet-Blondel method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.cc b/src/algorithms/reco/InclusiveKinematicsSigma.cc index 020d1dc5fb..c63b8f3981 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicsSigma.cc @@ -22,11 +22,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsSigma::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicsSigma::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -41,7 +40,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -55,7 +54,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -82,7 +81,7 @@ namespace eicrecon { auto gamma_h = hfs->at(0).getGamma(); if (sigma_h <= 0) { - m_log->debug("No scattered electron found or sigma zero or negative"); + debug("No scattered electron found or sigma zero or negative"); return; } @@ -97,7 +96,7 @@ namespace eicrecon { auto kin = kinematics->create(x_sig, Q2_sig, W_sig, y_sig, nu_sig); kin.setScat(kf); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsSigma.h b/src/algorithms/reco/InclusiveKinematicsSigma.h index d3849806c3..6cf3d13cbe 100644 --- a/src/algorithms/reco/InclusiveKinematicsSigma.h +++ b/src/algorithms/reco/InclusiveKinematicsSigma.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsSigmaAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsSigma - : public InclusiveKinematicsSigmaAlgorithm { - - public: - InclusiveKinematicsSigma(std::string_view name) - : InclusiveKinematicsSigmaAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using Sigma method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsSigmaAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicsSigma : public InclusiveKinematicsSigmaAlgorithm { + +public: + InclusiveKinematicsSigma(std::string_view name) + : InclusiveKinematicsSigmaAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using Sigma method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.cc b/src/algorithms/reco/InclusiveKinematicsTruth.cc index 54bc232f88..e48bbc0d77 100644 --- a/src/algorithms/reco/InclusiveKinematicsTruth.cc +++ b/src/algorithms/reco/InclusiveKinematicsTruth.cc @@ -19,9 +19,7 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicsTruth::init(std::shared_ptr& logger) { - m_log = logger; - + void InclusiveKinematicsTruth::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { // error() << "Unable to locate Particle Service. " @@ -47,7 +45,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const auto ei_p = ei_coll[0].getMomentum(); @@ -58,7 +56,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const auto pi_p = pi_coll[0].getMomentum(); @@ -73,7 +71,7 @@ namespace eicrecon { // the beam. const auto ef_coll = find_first_scattered_electron(mcparts); if (ef_coll.size() == 0) { - m_log->debug("No truth scattered electron found"); + debug("No truth scattered electron found"); return; } const auto ef_p = ef_coll[0].getMomentum(); @@ -91,7 +89,7 @@ namespace eicrecon { const auto W = sqrt(pi_mass*pi_mass + 2.*q_dot_pi - Q2); auto kin = kinematics->create(x, Q2, W, y, nu); - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicsTruth.h b/src/algorithms/reco/InclusiveKinematicsTruth.h index 8ec7a58cf7..83a3f44576 100644 --- a/src/algorithms/reco/InclusiveKinematicsTruth.h +++ b/src/algorithms/reco/InclusiveKinematicsTruth.h @@ -6,39 +6,30 @@ #include #include #include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicsTruthAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicsTruth - : public InclusiveKinematicsTruthAlgorithm { - - public: - InclusiveKinematicsTruth(std::string_view name) - : InclusiveKinematicsTruthAlgorithm{name, - {"MCParticles"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics from truth information."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicsTruthAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +class InclusiveKinematicsTruth : public InclusiveKinematicsTruthAlgorithm { + +public: + InclusiveKinematicsTruth(std::string_view name) + : InclusiveKinematicsTruthAlgorithm{ + name, + {"MCParticles"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics from truth information."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.cc b/src/algorithms/reco/InclusiveKinematicseSigma.cc index b898e41da9..56fb00eec1 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.cc +++ b/src/algorithms/reco/InclusiveKinematicseSigma.cc @@ -23,11 +23,10 @@ using ROOT::Math::PxPyPzEVector; namespace eicrecon { - void InclusiveKinematicseSigma::init(std::shared_ptr& logger) { - m_log = logger; + void InclusiveKinematicseSigma::init() { // m_pidSvc = service("ParticleSvc"); // if (!m_pidSvc) { - // m_log->debug("Unable to locate Particle Service. " + // debug("Unable to locate Particle Service. " // "Make sure you have ParticleSvc in the configuration."); // } } @@ -42,7 +41,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcparts); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector ei( @@ -56,7 +55,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcparts); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector pi( @@ -83,7 +82,7 @@ namespace eicrecon { auto gamma_h = hfs->at(0).getGamma(); if (sigma_h <= 0) { - m_log->debug("No scattered electron found or sigma zero or negative"); + debug("No scattered electron found or sigma zero or negative"); return; } @@ -106,7 +105,7 @@ namespace eicrecon { kin.setScat(kf); // Debugging output - m_log->debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), + debug("x,Q2,W,y,nu = {},{},{},{},{}", kin.getX(), kin.getQ2(), kin.getW(), kin.getY(), kin.getNu()); } diff --git a/src/algorithms/reco/InclusiveKinematicseSigma.h b/src/algorithms/reco/InclusiveKinematicseSigma.h index 9a505f0eae..030a37453d 100644 --- a/src/algorithms/reco/InclusiveKinematicseSigma.h +++ b/src/algorithms/reco/InclusiveKinematicseSigma.h @@ -4,45 +4,35 @@ #pragma once #include +#include #include #include #include -#include -#include -#include #include #include - namespace eicrecon { - using InclusiveKinematicseSigmaAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::ReconstructedParticleCollection, - edm4eic::HadronicFinalStateCollection - >, - algorithms::Output< - edm4eic::InclusiveKinematicsCollection - > - >; - - class InclusiveKinematicseSigma - : public InclusiveKinematicseSigmaAlgorithm { - - public: - InclusiveKinematicseSigma(std::string_view name) - : InclusiveKinematicseSigmaAlgorithm{name, - {"MCParticles", "scatteredElectron", "hadronicFinalState"}, - {"inclusiveKinematics"}, - "Determine inclusive kinematics using e-Sigma method."} {} - - void init(std::shared_ptr& logger); - void process(const Input&, const Output&) const final; - - private: - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - }; +using InclusiveKinematicseSigmaAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; + +class InclusiveKinematicseSigma : public InclusiveKinematicseSigmaAlgorithm { + +public: + InclusiveKinematicseSigma(std::string_view name) + : InclusiveKinematicseSigmaAlgorithm{ + name, + {"MCParticles", "scatteredElectron", "hadronicFinalState"}, + {"inclusiveKinematics"}, + "Determine inclusive kinematics using e-Sigma method."} {} + + void init() final; + void process(const Input&, const Output&) const final; + +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; +}; } // namespace eicrecon diff --git a/src/algorithms/reco/JetReconstruction.cc b/src/algorithms/reco/JetReconstruction.cc index ebe60f55aa..09dad9e089 100644 --- a/src/algorithms/reco/JetReconstruction.cc +++ b/src/algorithms/reco/JetReconstruction.cc @@ -6,9 +6,10 @@ // for error handling #include -#include // IWYU pragma: keep +#include // IWYU pragma: keep #include #include +#include #include // for fastjet objects #include @@ -24,139 +25,141 @@ using namespace fastjet; namespace eicrecon { - template - void JetReconstruction::init(std::shared_ptr logger) { - - m_log = logger; - m_log->trace("Initialized"); - - // if specified algorithm, recomb. scheme, or area type - // are not defined, then issue error and throw exception - try { - m_mapJetAlgo.at(m_cfg.jetAlgo); - } catch (std::out_of_range &out) { - m_log->error(" Unknown jet algorithm \"{}\" specified!", m_cfg.jetAlgo); - throw JException(out.what()); - } - - try { - m_mapRecombScheme.at(m_cfg.recombScheme); - } catch (std::out_of_range &out) { - m_log->error(" Unknown recombination scheme \"{}\" specified!", m_cfg.recombScheme); - throw JException(out.what()); - } - - try { - m_mapAreaType.at(m_cfg.areaType); - } catch (std::out_of_range &out) { - m_log->error(" Unknown area type \"{}\" specified!", m_cfg.areaType); - throw JException(out.what()); +template void JetReconstruction::init() { + + this->trace("Initialized"); + + // if specified algorithm, recomb. scheme, or area type + // are not defined, then issue error and throw exception + try { + m_mapJetAlgo.at(m_cfg.jetAlgo); + } catch (std::out_of_range& out) { + this->error(" Unknown jet algorithm \"{}\" specified!", m_cfg.jetAlgo); + throw JException(out.what()); + } + + try { + m_mapRecombScheme.at(m_cfg.recombScheme); + } catch (std::out_of_range& out) { + this->error(" Unknown recombination scheme \"{}\" specified!", m_cfg.recombScheme); + throw JException(out.what()); + } + + try { + m_mapAreaType.at(m_cfg.areaType); + } catch (std::out_of_range& out) { + this->error(" Unknown area type \"{}\" specified!", m_cfg.areaType); + throw JException(out.what()); + } + + // Choose jet definition based on no. of parameters + switch (m_mapJetAlgo[m_cfg.jetAlgo]) { + + // contributed algorithms + case JetAlgorithm::plugin_algorithm: + + // expand to other algorithms as required + if (m_cfg.jetContribAlgo == "Centauro") { + m_jet_plugin = std::make_unique(m_cfg.rJet); + m_jet_def = std::make_unique(m_jet_plugin.get()); + } else { + this->error(" Unknown contributed FastJet algorithm \"{}\" specified!", m_cfg.jetContribAlgo); + throw JException("Invalid contributed FastJet algorithm"); } - - // Choose jet definition based on no. of parameters - switch (m_mapJetAlgo[m_cfg.jetAlgo]) { - - // contributed algorithms - case JetAlgorithm::plugin_algorithm: - - // expand to other algorithms as required - if(m_cfg.jetContribAlgo == "Centauro"){ - m_jet_plugin = std::make_unique(m_cfg.rJet); - m_jet_def = std::make_unique(m_jet_plugin.get()); - } - else { - m_log->error(" Unknown contributed FastJet algorithm \"{}\" specified!", m_cfg.jetContribAlgo); - throw JException("Invalid contributed FastJet algorithm"); - } - break; - - // 0 parameter algorithms - case JetAlgorithm::ee_kt_algorithm: - m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_mapRecombScheme[m_cfg.recombScheme]); - break; - - // 2 parameter algorithms - case JetAlgorithm::genkt_algorithm: - [[fallthrough]]; - - case JetAlgorithm::ee_genkt_algorithm: - m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, m_cfg.pJet, m_mapRecombScheme[m_cfg.recombScheme]); - break; - - // all others have only 1 parameter - default: - m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, m_mapRecombScheme[m_cfg.recombScheme]); - break; - - } // end switch (jet algorithm) - - // Define jet area - m_area_def = std::make_unique(m_mapAreaType[m_cfg.areaType], GhostedAreaSpec(m_cfg.ghostMaxRap, m_cfg.numGhostRepeat, m_cfg.ghostArea)); - - } // end 'init(std::shared_ptr)' - - template - void JetReconstruction::process( - const typename JetReconstructionAlgorithm::Input& input, - const typename JetReconstructionAlgorithm::Output& output - ) const { - // Grab input collections - const auto [input_collection] = input; - auto [jet_collection] = output; - - // extract input momenta and collect into pseudojets - std::vector particles; - for (unsigned iInput = 0; const auto& input : *input_collection) { - - // get 4-vector - const auto& momentum = input.getMomentum(); - const auto& energy = input.getEnergy(); - const auto pt = edm4hep::utils::magnitudeTransverse(momentum); - - // Only cluster particles within the given pt Range - if ((pt > m_cfg.minCstPt) && (pt < m_cfg.maxCstPt)) { - particles.emplace_back(momentum.x, momentum.y, momentum.z, energy); - particles.back().set_user_index(iInput); - } - ++iInput; + break; + + // 0 parameter algorithms + case JetAlgorithm::ee_kt_algorithm: + m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], + m_mapRecombScheme[m_cfg.recombScheme]); + break; + + // 2 parameter algorithms + case JetAlgorithm::genkt_algorithm: + [[fallthrough]]; + + case JetAlgorithm::ee_genkt_algorithm: + m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, m_cfg.pJet, + m_mapRecombScheme[m_cfg.recombScheme]); + break; + + // all others have only 1 parameter + default: + m_jet_def = std::make_unique(m_mapJetAlgo[m_cfg.jetAlgo], m_cfg.rJet, + m_mapRecombScheme[m_cfg.recombScheme]); + break; + + } // end switch (jet algorithm) + + // Define jet area + m_area_def = std::make_unique( + m_mapAreaType[m_cfg.areaType], + GhostedAreaSpec(m_cfg.ghostMaxRap, m_cfg.numGhostRepeat, m_cfg.ghostArea)); + +} // end 'init()' + +template +void JetReconstruction::process( + const typename JetReconstructionAlgorithm::Input& input, + const typename JetReconstructionAlgorithm::Output& output) const { + // Grab input collections + const auto [input_collection] = input; + auto [jet_collection] = output; + + // extract input momenta and collect into pseudojets + std::vector particles; + for (unsigned iInput = 0; const auto& input : *input_collection) { + + // get 4-vector + const auto& momentum = input.getMomentum(); + const auto& energy = input.getEnergy(); + const auto pt = edm4hep::utils::magnitudeTransverse(momentum); + + // Only cluster particles within the given pt Range + if ((pt > m_cfg.minCstPt) && (pt < m_cfg.maxCstPt)) { + particles.emplace_back(momentum.x, momentum.y, momentum.z, energy); + particles.back().set_user_index(iInput); } + ++iInput; + } - // Skip empty - if (particles.empty()) { - m_log->trace(" Empty particle list."); - return; - } - m_log->trace(" Number of particles: {}", particles.size()); + // Skip empty + if (particles.empty()) { + this->trace(" Empty particle list."); + return; + } + this->trace(" Number of particles: {}", particles.size()); - // Run the clustering, extract the jets - fastjet::ClusterSequenceArea m_clus_seq(particles, *m_jet_def, *m_area_def); - std::vector jets = sorted_by_pt(m_clus_seq.inclusive_jets(m_cfg.minJetPt)); + // Run the clustering, extract the jets + fastjet::ClusterSequenceArea m_clus_seq(particles, *m_jet_def, *m_area_def); + std::vector jets = sorted_by_pt(m_clus_seq.inclusive_jets(m_cfg.minJetPt)); - // Print out some infos - m_log->trace(" Clustering with : {}", m_jet_def->description()); + // Print out some infos + this->trace(" Clustering with : {}", m_jet_def->description()); - // loop over jets - for (unsigned i = 0; i < jets.size(); i++) { + // loop over jets + for (unsigned i = 0; i < jets.size(); i++) { - m_log->trace(" jet {}: pt = {}, y = {}, phi = {}", i, jets[i].pt(), jets[i].rap(), jets[i].phi()); + this->trace(" jet {}: pt = {}, y = {}, phi = {}", i, jets[i].pt(), jets[i].rap(), + jets[i].phi()); - // create jet to store in output collection - edm4eic::MutableReconstructedParticle jet_output = jet_collection->create(); - jet_output.setMomentum(edm4hep::Vector3f(jets[i].px(), jets[i].py(), jets[i].pz())); - jet_output.setEnergy(jets[i].e()); - jet_output.setMass(jets[i].m()); + // create jet to store in output collection + edm4eic::MutableReconstructedParticle jet_output = jet_collection->create(); + jet_output.setMomentum(edm4hep::Vector3f(jets[i].px(), jets[i].py(), jets[i].pz())); + jet_output.setEnergy(jets[i].e()); + jet_output.setMass(jets[i].m()); - // link constituents to jet kinematic info - std::vector csts = jets[i].constituents(); - for (unsigned j = 0; j < csts.size(); j++) { - jet_output.addToParticles(input_collection->at(csts[j].user_index())); - } // for constituent j - } // for jet i + // link constituents to jet kinematic info + std::vector csts = jets[i].constituents(); + for (unsigned j = 0; j < csts.size(); j++) { + jet_output.addToParticles(input_collection->at(csts[j].user_index())); + } // for constituent j + } // for jet i - // return the jets - return; - } // end 'process(const T&)' + // return the jets + return; +} // end 'process(const T&)' - template class JetReconstruction; +template class JetReconstruction; -} // end namespace eicrecon +} // end namespace eicrecon diff --git a/src/algorithms/reco/JetReconstruction.h b/src/algorithms/reco/JetReconstruction.h index 460aa1d33d..9b2182479d 100644 --- a/src/algorithms/reco/JetReconstruction.h +++ b/src/algorithms/reco/JetReconstruction.h @@ -6,9 +6,7 @@ #include #include #include -#include #include -#include #include #include #include @@ -20,87 +18,73 @@ namespace eicrecon { - template - using JetReconstructionAlgorithm = algorithms::Algorithm< - algorithms::Input< - typename InputT::collection_type - >, - algorithms::Output< - edm4eic::ReconstructedParticleCollection - > - >; - - template - class JetReconstruction - : public JetReconstructionAlgorithm, - public WithPodConfig { - - public: - - JetReconstruction(std::string_view name) : - JetReconstructionAlgorithm { - name, - {"inputReconstructedParticles"}, - {"outputReconstructedParticles"}, - "Performs jet reconstruction using a FastJet algorithm." - } {} - - public: - - // algorithm initialization - void init(std::shared_ptr logger); - - // run algorithm - void process(const typename eicrecon::JetReconstructionAlgorithm::Input&, const typename eicrecon::JetReconstructionAlgorithm::Output&) const final; - - private: - - std::shared_ptr m_log; - - // fastjet components - std::unique_ptr m_jet_def; - std::unique_ptr m_area_def; - std::unique_ptr m_jet_plugin; - - // maps of user input onto fastjet options - std::map m_mapJetAlgo = { - {"kt_algorithm", fastjet::JetAlgorithm::kt_algorithm}, - {"cambridge_algorithm", fastjet::JetAlgorithm::cambridge_algorithm}, - {"antikt_algorithm", fastjet::JetAlgorithm::antikt_algorithm}, - {"genkt_algorithm", fastjet::JetAlgorithm::genkt_algorithm}, - {"cambridge_for_passive_algorithm", fastjet::JetAlgorithm::cambridge_for_passive_algorithm}, - {"genkt_for_passive_algorithm", fastjet::JetAlgorithm::genkt_for_passive_algorithm}, - {"ee_kt_algorithm", fastjet::JetAlgorithm::ee_kt_algorithm}, - {"ee_genkt_algorithm", fastjet::JetAlgorithm::ee_genkt_algorithm}, - {"plugin_algorithm", fastjet::JetAlgorithm::plugin_algorithm} - }; - std::map m_mapRecombScheme = { - {"E_scheme", fastjet::RecombinationScheme::E_scheme}, - {"pt_scheme", fastjet::RecombinationScheme::pt_scheme}, - {"pt2_scheme", fastjet::RecombinationScheme::pt2_scheme}, - {"Et_scheme", fastjet::RecombinationScheme::Et_scheme}, - {"Et2_scheme", fastjet::RecombinationScheme::Et2_scheme}, - {"BIpt_scheme", fastjet::RecombinationScheme::BIpt_scheme}, - {"BIpt2_scheme", fastjet::RecombinationScheme::BIpt2_scheme}, - {"WTA_pt_scheme", fastjet::RecombinationScheme::WTA_pt_scheme}, - {"WTA_modp_scheme", fastjet::RecombinationScheme::WTA_modp_scheme}, - {"external_scheme", fastjet::RecombinationScheme::external_scheme} - }; - std::map m_mapAreaType = { - {"active_area", fastjet::AreaType::active_area}, - {"active_area_explicit_ghosts", fastjet::AreaType::active_area_explicit_ghosts}, - {"one_ghost_passive_area", fastjet::AreaType::one_ghost_passive_area}, - {"passive_area", fastjet::AreaType::passive_area}, - {"voronoi_area", fastjet::AreaType::voronoi_area} - }; - - // default fastjet options - const struct defaults { - std::string jetAlgo; - std::string recombScheme; - std::string areaType; - } m_defaultFastjetOpts = {"antikt_algorithm", "E_scheme", "active_area"}; - - }; // end JetReconstruction definition - -} // end eicrecon namespace +template +using JetReconstructionAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +template +class JetReconstruction : public JetReconstructionAlgorithm, + public WithPodConfig { + +public: + JetReconstruction(std::string_view name) + : JetReconstructionAlgorithm{ + name, + {"inputReconstructedParticles"}, + {"outputReconstructedParticles"}, + "Performs jet reconstruction using a FastJet algorithm."} {} + +public: + // algorithm initialization + void init() final; + + // run algorithm + void process(const typename eicrecon::JetReconstructionAlgorithm::Input&, + const typename eicrecon::JetReconstructionAlgorithm::Output&) const final; + +private: + // fastjet components + std::unique_ptr m_jet_def; + std::unique_ptr m_area_def; + std::unique_ptr m_jet_plugin; + + // maps of user input onto fastjet options + std::map m_mapJetAlgo = { + {"kt_algorithm", fastjet::JetAlgorithm::kt_algorithm}, + {"cambridge_algorithm", fastjet::JetAlgorithm::cambridge_algorithm}, + {"antikt_algorithm", fastjet::JetAlgorithm::antikt_algorithm}, + {"genkt_algorithm", fastjet::JetAlgorithm::genkt_algorithm}, + {"cambridge_for_passive_algorithm", fastjet::JetAlgorithm::cambridge_for_passive_algorithm}, + {"genkt_for_passive_algorithm", fastjet::JetAlgorithm::genkt_for_passive_algorithm}, + {"ee_kt_algorithm", fastjet::JetAlgorithm::ee_kt_algorithm}, + {"ee_genkt_algorithm", fastjet::JetAlgorithm::ee_genkt_algorithm}, + {"plugin_algorithm", fastjet::JetAlgorithm::plugin_algorithm}}; + std::map m_mapRecombScheme = { + {"E_scheme", fastjet::RecombinationScheme::E_scheme}, + {"pt_scheme", fastjet::RecombinationScheme::pt_scheme}, + {"pt2_scheme", fastjet::RecombinationScheme::pt2_scheme}, + {"Et_scheme", fastjet::RecombinationScheme::Et_scheme}, + {"Et2_scheme", fastjet::RecombinationScheme::Et2_scheme}, + {"BIpt_scheme", fastjet::RecombinationScheme::BIpt_scheme}, + {"BIpt2_scheme", fastjet::RecombinationScheme::BIpt2_scheme}, + {"WTA_pt_scheme", fastjet::RecombinationScheme::WTA_pt_scheme}, + {"WTA_modp_scheme", fastjet::RecombinationScheme::WTA_modp_scheme}, + {"external_scheme", fastjet::RecombinationScheme::external_scheme}}; + std::map m_mapAreaType = { + {"active_area", fastjet::AreaType::active_area}, + {"active_area_explicit_ghosts", fastjet::AreaType::active_area_explicit_ghosts}, + {"one_ghost_passive_area", fastjet::AreaType::one_ghost_passive_area}, + {"passive_area", fastjet::AreaType::passive_area}, + {"voronoi_area", fastjet::AreaType::voronoi_area}}; + + // default fastjet options + const struct defaults { + std::string jetAlgo; + std::string recombScheme; + std::string areaType; + } m_defaultFastjetOpts = {"antikt_algorithm", "E_scheme", "active_area"}; + +}; // end JetReconstruction definition + +} // namespace eicrecon diff --git a/src/algorithms/reco/TransformBreitFrame.cc b/src/algorithms/reco/TransformBreitFrame.cc index c8eff8221d..5e168845d3 100644 --- a/src/algorithms/reco/TransformBreitFrame.cc +++ b/src/algorithms/reco/TransformBreitFrame.cc @@ -20,14 +20,6 @@ namespace eicrecon { - void TransformBreitFrame::init(std::shared_ptr logger) { - - m_log = logger; - m_log->trace("Initialized"); - - } // end 'init(std::shared_ptr)' - - void TransformBreitFrame::process( const TransformBreitFrame::Input& input, const TransformBreitFrame::Output& output @@ -42,7 +34,7 @@ namespace eicrecon { // Get incoming electron beam const auto ei_coll = find_first_beam_electron(mcpart); if (ei_coll.size() == 0) { - m_log->debug("No beam electron found"); + debug("No beam electron found"); return; } const PxPyPzEVector e_initial( @@ -56,7 +48,7 @@ namespace eicrecon { // Get incoming hadron beam const auto pi_coll = find_first_beam_hadron(mcpart); if (pi_coll.size() == 0) { - m_log->debug("No beam hadron found"); + debug("No beam hadron found"); return; } const PxPyPzEVector p_initial( @@ -67,11 +59,11 @@ namespace eicrecon { m_crossingAngle) ); - m_log->debug("electron energy, proton energy = {},{}",e_initial.E(),p_initial.E()); + debug("electron energy, proton energy = {},{}",e_initial.E(),p_initial.E()); // Get the event kinematics, set up transform if (kine->size() == 0) { - m_log->debug("No kinematics found"); + debug("No kinematics found"); return; } @@ -82,15 +74,15 @@ namespace eicrecon { // Use relation to get reconstructed scattered electron const PxPyPzEVector e_final = edm4hep::utils::detail::p4(evt_kin.getScat(),&edm4hep::utils::UseEnergy); - m_log->debug("scattered electron in lab frame px,py,pz,E = {},{},{},{}", + debug("scattered electron in lab frame px,py,pz,E = {},{},{},{}", e_final.Px(),e_final.Py(),e_final.Pz(),e_final.E()); // Set up the transformation const PxPyPzEVector virtual_photon = (e_initial - e_final); - m_log->debug("virtual photon in lab frame px,py,pz,E = {},{},{},{}", + debug("virtual photon in lab frame px,py,pz,E = {},{},{},{}", virtual_photon.Px(),virtual_photon.Py(),virtual_photon.Pz(),virtual_photon.E()); - m_log->debug("x, Q^2 = {},{}",meas_x,meas_Q2); + debug("x, Q^2 = {},{}",meas_x,meas_Q2); // Set up the transformation (boost) to the Breit frame const auto P3 = p_initial.Vect(); @@ -116,9 +108,9 @@ namespace eicrecon { e_final_breit = breitRot*e_final_breit; virtual_photon_breit = breitRot*virtual_photon_breit; - m_log->debug("incoming hadron in Breit frame px,py,pz,E = {},{},{},{}", + debug("incoming hadron in Breit frame px,py,pz,E = {},{},{},{}", p_initial_breit.Px(),p_initial_breit.Py(),p_initial_breit.Pz(),p_initial_breit.E()); - m_log->debug("virtual photon in Breit frame px,py,pz,E = {},{},{},{}", + debug("virtual photon in Breit frame px,py,pz,E = {},{},{},{}", virtual_photon_breit.Px(),virtual_photon_breit.Py(),virtual_photon_breit.Pz(),virtual_photon_breit.E()); // look over the input particles and transform diff --git a/src/algorithms/reco/TransformBreitFrame.h b/src/algorithms/reco/TransformBreitFrame.h index af654c2772..74a89dd69b 100644 --- a/src/algorithms/reco/TransformBreitFrame.h +++ b/src/algorithms/reco/TransformBreitFrame.h @@ -7,8 +7,6 @@ #include #include #include -#include -#include #include #include @@ -16,42 +14,30 @@ namespace eicrecon { - using TransformBreitFrameAlgorithm = algorithms::Algorithm< - algorithms::Input< - edm4hep::MCParticleCollection, - edm4eic::InclusiveKinematicsCollection, - edm4eic::ReconstructedParticleCollection - >, - algorithms::Output< - edm4eic::ReconstructedParticleCollection - > - >; +using TransformBreitFrameAlgorithm = algorithms::Algorithm< + algorithms::Input, + algorithms::Output>; - class TransformBreitFrame - : public TransformBreitFrameAlgorithm, - public WithPodConfig { +class TransformBreitFrame : public TransformBreitFrameAlgorithm, public WithPodConfig { - public: +public: + TransformBreitFrame(std::string_view name) + : TransformBreitFrameAlgorithm{ + name, + {"inputMCParticles", "inputInclusiveKinematics", "inputReconstructedParticles"}, + {"outputReconstructedParticles"}, + "Transforms a set of particles from the lab frame to the Breit frame"} {} - TransformBreitFrame(std::string_view name) : - TransformBreitFrameAlgorithm { - name, - {"inputMCParticles", "inputInclusiveKinematics", "inputReconstructedParticles"}, - {"outputReconstructedParticles"}, - "Transforms a set of particles from the lab frame to the Breit frame" - } {} + // algorithm initialization + void init() final{}; - // algorithm initialization - void init(std::shared_ptr logger); + // run algorithm + void process(const Input&, const Output&) const final; - // run algorithm - void process(const Input&, const Output&) const final; +private: + double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - private: +}; // end TransformBreitFrame definition - std::shared_ptr m_log; - double m_proton{0.93827}, m_neutron{0.93957}, m_electron{0.000510998928}, m_crossingAngle{-0.025}; - - }; // end TransformBreitFrame definition - -} // end eicrecon namespace +} // namespace eicrecon diff --git a/src/factories/digi/SiliconTrackerDigi_factory.h b/src/factories/digi/SiliconTrackerDigi_factory.h index 3975dd34da..a708b3335c 100644 --- a/src/factories/digi/SiliconTrackerDigi_factory.h +++ b/src/factories/digi/SiliconTrackerDigi_factory.h @@ -30,7 +30,7 @@ class SiliconTrackerDigi_factory : public JOmniFactory(GetPrefix()); m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/fardetectors/FarDetectorLinearProjection_factory.h b/src/factories/fardetectors/FarDetectorLinearProjection_factory.h index 39510d37ca..ebaeab6539 100644 --- a/src/factories/fardetectors/FarDetectorLinearProjection_factory.h +++ b/src/factories/fardetectors/FarDetectorLinearProjection_factory.h @@ -29,7 +29,7 @@ public JOmniFactory(GetPrefix()); - m_algo->level((algorithms::LogLevel)logger()->level()); + m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); m_algo->init(); } diff --git a/src/factories/fardetectors/FarDetectorLinearTracking_factory.h b/src/factories/fardetectors/FarDetectorLinearTracking_factory.h index aa3dd32247..5611e1fc26 100644 --- a/src/factories/fardetectors/FarDetectorLinearTracking_factory.h +++ b/src/factories/fardetectors/FarDetectorLinearTracking_factory.h @@ -33,7 +33,7 @@ class FarDetectorLinearTracking_factory : m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h b/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h index 13ace79e81..2134c064b2 100644 --- a/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h +++ b/src/factories/fardetectors/FarDetectorMLReconstruction_factory.h @@ -39,7 +39,7 @@ class FarDetectorMLReconstruction_factory : public: void Configure() { m_algo = std::make_unique(GetPrefix()); - m_algo->level((algorithms::LogLevel)logger()->level()); + m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); m_algo->init(); } diff --git a/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h b/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h index d45aa97e3e..fa4b400d63 100644 --- a/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h +++ b/src/factories/fardetectors/FarDetectorTrackerCluster_factory.h @@ -36,7 +36,7 @@ public JOmniFactorylevel(static_cast(logger()->level())); // Setup algorithm m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } diff --git a/src/factories/meta/CollectionCollector_factory.h b/src/factories/meta/CollectionCollector_factory.h index 9edf056f58..df5f6da39c 100644 --- a/src/factories/meta/CollectionCollector_factory.h +++ b/src/factories/meta/CollectionCollector_factory.h @@ -22,7 +22,7 @@ namespace eicrecon { void Configure() { m_algo = std::make_unique(this->GetPrefix()); m_algo->level(static_cast(this->logger()->level())); - m_algo->init(this->logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/HadronicFinalState_factory.h b/src/factories/reco/HadronicFinalState_factory.h index 30b0dea850..00d03523e3 100644 --- a/src/factories/reco/HadronicFinalState_factory.h +++ b/src/factories/reco/HadronicFinalState_factory.h @@ -32,7 +32,8 @@ class HadronicFinalState_factory : public: void Configure() { m_algo = std::make_unique(this->GetPrefix()); - m_algo->init(this->logger()); + m_algo->level(static_cast(this->logger()->level())); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/InclusiveKinematicsML_factory.h b/src/factories/reco/InclusiveKinematicsML_factory.h index 9f4ccd10d2..78343cfbd9 100644 --- a/src/factories/reco/InclusiveKinematicsML_factory.h +++ b/src/factories/reco/InclusiveKinematicsML_factory.h @@ -34,7 +34,7 @@ class InclusiveKinematicsML_factory : m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); m_algo->applyConfig(config()); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h index 7e99ba0253..7c395f39dc 100644 --- a/src/factories/reco/InclusiveKinematicsReconstructed_factory.h +++ b/src/factories/reco/InclusiveKinematicsReconstructed_factory.h @@ -33,7 +33,7 @@ class InclusiveKinematicsReconstructed_factory : void Configure() { m_algo = std::make_unique(this->GetPrefix()); m_algo->level(static_cast(this->logger()->level())); - m_algo->init(this->logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/InclusiveKinematicsTruth_factory.h b/src/factories/reco/InclusiveKinematicsTruth_factory.h index 6c5f644fdc..b119c0541b 100644 --- a/src/factories/reco/InclusiveKinematicsTruth_factory.h +++ b/src/factories/reco/InclusiveKinematicsTruth_factory.h @@ -30,7 +30,7 @@ class InclusiveKinematicsTruth_factory : void Configure() { m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/JetReconstruction_factory.h b/src/factories/reco/JetReconstruction_factory.h index 98a665e826..6919ae5d08 100644 --- a/src/factories/reco/JetReconstruction_factory.h +++ b/src/factories/reco/JetReconstruction_factory.h @@ -46,7 +46,7 @@ namespace eicrecon { m_algo = std::make_unique(this->GetPrefix()); m_algo->level(static_cast(this->logger()->level())); m_algo->applyConfig(FactoryT::config()); - m_algo->init(FactoryT::logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/factories/reco/TransformBreitFrame_factory.h b/src/factories/reco/TransformBreitFrame_factory.h index 70836b075b..3a9852818b 100644 --- a/src/factories/reco/TransformBreitFrame_factory.h +++ b/src/factories/reco/TransformBreitFrame_factory.h @@ -38,7 +38,7 @@ namespace eicrecon { void Configure() { m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { diff --git a/src/global/pid/MergeTrack_factory.h b/src/global/pid/MergeTrack_factory.h index e3fa78b518..c59f4250ba 100644 --- a/src/global/pid/MergeTrack_factory.h +++ b/src/global/pid/MergeTrack_factory.h @@ -33,7 +33,7 @@ class MergeTrack_factory : public JOmniFactory { void Configure() { m_algo = std::make_unique(GetPrefix()); m_algo->level(static_cast(logger()->level())); - m_algo->init(logger()); + m_algo->init(); } void ChangeRun(int64_t run_number) { } diff --git a/src/tests/algorithms_test/pid_MergeTracks.cc b/src/tests/algorithms_test/pid_MergeTracks.cc index 1bf10a8c09..ee2d3fd3e1 100644 --- a/src/tests/algorithms_test/pid_MergeTracks.cc +++ b/src/tests/algorithms_test/pid_MergeTracks.cc @@ -1,16 +1,15 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) 2023, Christopher Dilks +#include #include #include #include +#include #include #include #include -#include -#include -#include -#include +#include #include #include @@ -21,9 +20,8 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { // initialize algorithm //---------------------------------------------------------- - std::shared_ptr logger = spdlog::default_logger()->clone("MergeTracks"); - logger->set_level(spdlog::level::debug); - algo.init(logger); + algo.level(algorithms::LogLevel::kDebug); + algo.init(); // helper functions and objects //---------------------------------------------------------- @@ -32,12 +30,12 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { struct Point { decltype(edm4eic::TrackPoint::position) position; decltype(edm4eic::TrackPoint::momentum) momentum; - decltype(edm4eic::TrackPoint::time) time; + decltype(edm4eic::TrackPoint::time) time; }; - auto make_track = [] (auto& coll, std::vector points) { + auto make_track = [](auto& coll, std::vector points) { auto trk = coll->create(); - for(auto& point : points) { + for (auto& point : points) { edm4eic::TrackPoint trk_point; trk_point.position = point.position; trk_point.momentum = point.momentum; @@ -46,14 +44,11 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { } }; - auto track_length = [] (auto trk) { + auto track_length = [](auto trk) { auto a = trk.getPoints(0); - auto b = trk.getPoints(trk.points_size()-1); - return std::hypot( - a.position.x - b.position.x, - a.position.y - b.position.y, - a.position.z - b.position.z - ); + auto b = trk.getPoints(trk.points_size() - 1); + return std::hypot(a.position.x - b.position.x, a.position.y - b.position.y, + a.position.z - b.position.z); }; // inputs @@ -73,112 +68,101 @@ TEST_CASE("the PID MergeTracks algorithm runs", "[MergeTracks]") { auto collection2 = std::make_unique(); // track0 - make_track(collection0, { // horizontal - { {0, 0, 0}, {1, 0, 0}, 1 }, // { position, momentum, time } - { {1, 0, 0}, {1, 0, 0}, 2 } - }); - make_track(collection1, { // horizontal - { {2, 0, 0}, {1, 0, 0}, 3 }, - { {3, 0, 0}, {1, 0, 0}, 4 } - }); - make_track(collection2, { // horizontal - { {4, 0, 0}, {1, 0, 0}, 5 }, - { {5, 0, 0}, {1, 0, 0}, 6 } - }); + make_track(collection0, { // horizontal + {{0, 0, 0}, {1, 0, 0}, 1}, // { position, momentum, time } + {{1, 0, 0}, {1, 0, 0}, 2}}); + make_track(collection1, {// horizontal + {{2, 0, 0}, {1, 0, 0}, 3}, + {{3, 0, 0}, {1, 0, 0}, 4}}); + make_track(collection2, {// horizontal + {{4, 0, 0}, {1, 0, 0}, 5}, + {{5, 0, 0}, {1, 0, 0}, 6}}); // track1 - make_track(collection0, { // horizontal - { {0, 0, 0}, {1, 0, 0}, 1 }, - { {1, 0, 0}, {1, 0, 0}, 2 } - }); - make_track(collection1, { // empty - { } - }); - make_track(collection2, { // one point - { {2, 0, 0}, {1, 0, 0}, 3 } - }); + make_track(collection0, {// horizontal + {{0, 0, 0}, {1, 0, 0}, 1}, + {{1, 0, 0}, {1, 0, 0}, 2}}); + make_track(collection1, {// empty + {}}); + make_track(collection2, {// one point + {{2, 0, 0}, {1, 0, 0}, 3}}); // track2 - make_track(collection0, { // vertical - { {0, 0, 0}, {0, 1, 0}, 1 }, - { {0, 1, 0}, {0, 1, 0}, 2 } - }); - make_track(collection1, { // horizontal - { {0, 1, 0}, {1, 0, 0}, 3 }, - { {1, 1, 0}, {1, 0, 0}, 4 } - }); - make_track(collection2, { // vertical - { {1, 1, 0}, {0, 1, 0}, 5 }, - { {1, 2, 0}, {0, 1, 0}, 6 } - }); + make_track(collection0, {// vertical + {{0, 0, 0}, {0, 1, 0}, 1}, + {{0, 1, 0}, {0, 1, 0}, 2}}); + make_track(collection1, {// horizontal + {{0, 1, 0}, {1, 0, 0}, 3}, + {{1, 1, 0}, {1, 0, 0}, 4}}); + make_track(collection2, {// vertical + {{1, 1, 0}, {0, 1, 0}, 5}, + {{1, 2, 0}, {0, 1, 0}, 6}}); // track3 - make_track(collection0, { // horizontal - { {1, 0, 0}, {1, 0, 0}, 2 }, - { {0, 0, 0}, {1, 0, 0}, 1 } - }); - make_track(collection1, { // horizontal - { {3, 0, 0}, {1, 0, 0}, 4 }, - { {4, 0, 0}, {1, 0, 0}, 5 } - }); - make_track(collection2, { // horizontal - { {5, 0, 0}, {1, 0, 0}, 6 }, - { {2, 0, 0}, {1, 0, 0}, 3 } - }); + make_track(collection0, {// horizontal + {{1, 0, 0}, {1, 0, 0}, 2}, + {{0, 0, 0}, {1, 0, 0}, 1}}); + make_track(collection1, {// horizontal + {{3, 0, 0}, {1, 0, 0}, 4}, + {{4, 0, 0}, {1, 0, 0}, 5}}); + make_track(collection2, {// horizontal + {{5, 0, 0}, {1, 0, 0}, 6}, + {{2, 0, 0}, {1, 0, 0}, 3}}); // tests //---------------------------------------------------------- SECTION("merge tracks from 1 collection") { // run algorithm - std::vector> colls = { collection0.get() }; + std::vector> colls = {collection0.get()}; // create output collection auto trks = std::make_unique(); algo.process({colls}, {trks.get()}); // input collection(s) size == output collection size REQUIRE(trks->size() == colls.front()->size()); // track length: from endpoints - REQUIRE_THAT( track_length(trks->at(0)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(2)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(3)), Catch::Matchers::WithinAbs(1, EPSILON) ); + REQUIRE_THAT(track_length(trks->at(0)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(2)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(3)), Catch::Matchers::WithinAbs(1, EPSILON)); // track length: from algorithm // FIXME when implemented in `MergeTracks` - for(const auto& trk : *trks) - REQUIRE_THAT( trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON) ); + for (const auto& trk : *trks) + REQUIRE_THAT(trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON)); } SECTION("merge tracks from 2 collections") { // run algorithm - std::vector> colls = { collection0.get(), collection1.get() }; + std::vector> colls = {collection0.get(), + collection1.get()}; auto trks = std::make_unique(); algo.process({colls}, {trks.get()}); // input collection(s) size == output collection size REQUIRE(trks->size() == colls.front()->size()); // track length: from endpoints - REQUIRE_THAT( track_length(trks->at(0)), Catch::Matchers::WithinAbs(3, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1,1), EPSILON) ); - REQUIRE_THAT( track_length(trks->at(3)), Catch::Matchers::WithinAbs(4, EPSILON) ); + REQUIRE_THAT(track_length(trks->at(0)), Catch::Matchers::WithinAbs(3, EPSILON)); + REQUIRE_THAT(track_length(trks->at(1)), Catch::Matchers::WithinAbs(1, EPSILON)); + REQUIRE_THAT(track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1, 1), EPSILON)); + REQUIRE_THAT(track_length(trks->at(3)), Catch::Matchers::WithinAbs(4, EPSILON)); // track length: from algorithm // FIXME when implemented in `MergeTracks` - for(const auto& trk : *trks) - REQUIRE_THAT( trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON) ); + for (const auto& trk : *trks) + REQUIRE_THAT(trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON)); } SECTION("merge tracks from 3 collections") { // run algorithm - std::vector> colls = { collection0.get(), collection1.get(), collection2.get() }; + std::vector> colls = { + collection0.get(), collection1.get(), collection2.get()}; auto trks = std::make_unique(); algo.process({colls}, {trks.get()}); // input collection(s) size == output collection size REQUIRE(trks->size() == colls.front()->size()); // track length: from endpoints - REQUIRE_THAT( track_length(trks->at(0)), Catch::Matchers::WithinAbs(5, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(1)), Catch::Matchers::WithinAbs(2, EPSILON) ); - REQUIRE_THAT( track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1,2), EPSILON) ); - REQUIRE_THAT( track_length(trks->at(3)), Catch::Matchers::WithinAbs(5, EPSILON) ); + REQUIRE_THAT(track_length(trks->at(0)), Catch::Matchers::WithinAbs(5, EPSILON)); + REQUIRE_THAT(track_length(trks->at(1)), Catch::Matchers::WithinAbs(2, EPSILON)); + REQUIRE_THAT(track_length(trks->at(2)), Catch::Matchers::WithinAbs(std::hypot(1, 2), EPSILON)); + REQUIRE_THAT(track_length(trks->at(3)), Catch::Matchers::WithinAbs(5, EPSILON)); // track length: from algorithm // FIXME when implemented in `MergeTracks` - for(const auto& trk : *trks) - REQUIRE_THAT( trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON) ); + for (const auto& trk : *trks) + REQUIRE_THAT(trk.getLength(), Catch::Matchers::WithinAbs(0, EPSILON)); } - } diff --git a/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc b/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc index 64c612782c..fd66879028 100644 --- a/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc +++ b/src/tests/algorithms_test/tracking_SiliconSimpleCluster.cc @@ -5,160 +5,143 @@ #include #include #include +#include #include #include #include #include -#include -#include -#include -#include #include -#include +#include #include #include #include "algorithms/fardetectors/FarDetectorTrackerCluster.h" #include "algorithms/fardetectors/FarDetectorTrackerClusterConfig.h" -TEST_CASE( "the clustering algorithm runs", "[FarDetectorTrackerCluster]" ) { +TEST_CASE("the clustering algorithm runs", "[FarDetectorTrackerCluster]") { eicrecon::FarDetectorTrackerCluster algo("FarDetectorTrackerCluster"); - std::shared_ptr logger = spdlog::default_logger()->clone("FarDetectorTrackerCluster"); - logger->set_level(spdlog::level::trace); - eicrecon::FarDetectorTrackerClusterConfig cfg; cfg.hit_time_limit = 10.0 * edm4eic::unit::ns; - cfg.readout = "MockTrackerHits"; - cfg.x_field = "x"; - cfg.y_field = "y"; + cfg.readout = "MockTrackerHits"; + cfg.x_field = "x"; + cfg.y_field = "y"; auto detector = algorithms::GeoSvc::instance().detector(); - auto id_desc = detector->readout(cfg.readout).idSpec(); + auto id_desc = detector->readout(cfg.readout).idSpec(); algo.applyConfig(cfg); - algo.init(logger); + algo.level(algorithms::LogLevel::kTrace); + algo.init(); - SECTION( "on a single pixel" ) { + SECTION("on a single pixel") { edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 0.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 0.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 1 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[0].x == 0.0 ); - REQUIRE( clusterPositions[0].y == 0.0 ); - REQUIRE( clusterPositions[0].energy == 5.0 ); - REQUIRE( clusterPositions[0].time == 0.0 ); + REQUIRE(clusterPositions.size() == 1); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[0].x == 0.0); + REQUIRE(clusterPositions[0].y == 0.0); + REQUIRE(clusterPositions[0].energy == 5.0); + REQUIRE(clusterPositions[0].time == 0.0); } - SECTION( "on two separated pixels" ) { + SECTION("on two separated pixels") { edm4eic::RawTrackerHitCollection hits_coll; hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 10}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + id_desc.encode({{"system", 255}, {"x", 0}, {"y", 10}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 10}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + id_desc.encode({{"system", 255}, {"x", 10}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 2 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[1].rawHits.size() == 1 ); + REQUIRE(clusterPositions.size() == 2); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[1].rawHits.size() == 1); } - SECTION( "on two adjacent pixels" ) { + SECTION("on two adjacent pixels") { edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 5.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 5.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 1 ); - REQUIRE( clusterPositions[0].rawHits.size() == 2 ); - REQUIRE( clusterPositions[0].x == 0.5 ); - + REQUIRE(clusterPositions.size() == 1); + REQUIRE(clusterPositions[0].rawHits.size() == 2); + REQUIRE(clusterPositions[0].x == 0.5); } - SECTION( "on two adjacent pixels outwith the time separation" ) { + SECTION("on two adjacent pixels outwith the time separation") { edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 0.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 0.0 // int32 timeStamp ); - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, - 5.0, // int32 charge, - 1.1 * cfg.hit_time_limit // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, + 5.0, // int32 charge, + 1.1 * cfg.hit_time_limit // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - REQUIRE( clusterPositions.size() == 2 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[1].rawHits.size() == 1 ); + REQUIRE(clusterPositions.size() == 2); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[1].rawHits.size() == 1); } - - SECTION( "run on three adjacent pixels" ) { + SECTION("run on three adjacent pixels") { // Check I and L shape clusters - auto pixel3 = GENERATE(std::vector{2,0}, std::vector{1,1}); - auto pixelCharges = GENERATE(std::vector{5,10,5}, std::vector{10,5,5}); - float pixel2Time = GENERATE_COPY(0, 1.1 * cfg.hit_time_limit); + auto pixel3 = GENERATE(std::vector{2, 0}, std::vector{1, 1}); + auto pixelCharges = GENERATE(std::vector{5, 10, 5}, std::vector{10, 5, 5}); + float pixel2Time = GENERATE_COPY(0, 1.1 * cfg.hit_time_limit); edm4eic::RawTrackerHitCollection hits_coll; - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, - pixelCharges[0], // int32 charge, - 0.0 // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}}), // std::uint64_t cellID, + pixelCharges[0], // int32 charge, + 0.0 // int32 timeStamp ); - hits_coll.create( - id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, - pixelCharges[1], // int32 charge, - pixel2Time // int32 timeStamp + hits_coll.create(id_desc.encode({{"system", 255}, {"x", 1}, {"y", 0}}), // std::uint64_t cellID, + pixelCharges[1], // int32 charge, + pixel2Time // int32 timeStamp ); hits_coll.create( - id_desc.encode({{"system", 255}, {"x", pixel3[0]}, {"y", pixel3[1]}}), // std::uint64_t cellID, - pixelCharges[2], // int32 charge, - 0.0 // int32 timeStamp + id_desc.encode( + {{"system", 255}, {"x", pixel3[0]}, {"y", pixel3[1]}}), // std::uint64_t cellID, + pixelCharges[2], // int32 charge, + 0.0 // int32 timeStamp ); std::vector clusterPositions = algo.ClusterHits(hits_coll); - if(pixel2Time < cfg.hit_time_limit) { - REQUIRE( clusterPositions.size() == 1 ); - REQUIRE( clusterPositions[0].rawHits.size() == 3 ); + if (pixel2Time < cfg.hit_time_limit) { + REQUIRE(clusterPositions.size() == 1); + REQUIRE(clusterPositions[0].rawHits.size() == 3); + } else if (pixel3[0] == 2) { + REQUIRE(clusterPositions.size() == 3); + REQUIRE(clusterPositions[0].rawHits.size() == 1); + REQUIRE(clusterPositions[1].rawHits.size() == 1); + REQUIRE(clusterPositions[2].rawHits.size() == 1); + } else { + REQUIRE(clusterPositions.size() == 2); } - else if(pixel3[0] == 2){ - REQUIRE( clusterPositions.size() == 3 ); - REQUIRE( clusterPositions[0].rawHits.size() == 1 ); - REQUIRE( clusterPositions[1].rawHits.size() == 1 ); - REQUIRE( clusterPositions[2].rawHits.size() == 1 ); - } - else { - REQUIRE( clusterPositions.size() == 2 ); - } - - } }