Commit 661a9c7c authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Added simple clustering algorithm

- simple clustering algorithm returns one cluster.
- uses 3D hits
parent 5f6e92d1
......@@ -33,6 +33,7 @@ gaudi_add_module(JugRecoPlugins
src/components/TrackFindingAlgorithm.cpp
src/components/TestACTSLogger.cpp
src/components/TrackParamTruthInit.cpp
src/components/SimpleClustering.cpp
src/components/ParticlesFromTrackFit.cpp
src/components/CalorimeterIslandCluster.cpp
src/components/CrystalEndcapsReco.cpp
......
......@@ -71,7 +71,7 @@ namespace Jug {
return;
}
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / params[Acts::eBoundQOverP], 0.000511},
eic::Particle p({params[Acts::eBoundPhi], params[Acts::eBoundTheta], 1.0 / std::abs(params[Acts::eBoundQOverP]), 0.000511},
{0.0, 0.0, 0.0, params[Acts::eBoundTime]},
(long long)11 * params[Acts::eBoundQOverP] / std::abs(params[Acts::eBoundQOverP]), 0);
//debug() << p << endmsg;
......
#include <algorithm>
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
// FCCSW
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
#include "eicd/ClusterCollection.h"
using namespace Gaudi::Units;
namespace Jug::Reco {
class SimpleClustering : public GaudiAlgorithm {
public:
using RecHits = eic::CalorimeterHitCollection;
using Clusters = eic::ClusterCollection;
DataHandle<RecHits> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<Clusters> m_outputClusters{"outputClusters", Gaudi::DataHandle::Writer, this};
Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 5.0*MeV};
Gaudi::Property<double> m_maxDistance{this, "maxDistance", 20.0*cm};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
SimpleClustering(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputClusters", m_outputClusters, "Output clusters");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
m_geoSvc = service("GeoSvc");
if (!m_geoSvc) {
error() << "Unable to locate Geometry Service. "
<< "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collections
const auto& hits = *m_inputHitCollection.get();
// Create output collections
auto& clusters = *m_outputClusters.createAndPut();
std::vector<eic::CalorimeterHit> first_cluster_hits ;
std::vector<eic::CalorimeterHit> remaining_hits ;
double max_dist = m_maxDistance.value();
debug() << " max_dist = " << max_dist << endmsg;
eic::CalorimeterHit ref_hit;
ref_hit.energy(0.0);
for (const auto& h : hits) {
if(h.energy() > ref_hit.energy()){
ref_hit = h;
}
}
debug() << " ref_hit" << ref_hit << endmsg;
//auto ref_hit = *(std::max_element(std::begin(hits), std::end(hits),
// [](const auto& h1, const auto& h2) { return h1.energy() < h2.energy(); }));
//std::partition_copy(std::begin(hits), std::end(hits), std::begin(first_cluster_hits), std::begin(remaining_hits),
// [=](const auto& h) {
// return (std::hypot(h.x() - ref_hit.x(), h.y() - ref_hit.y(), h.z() - ref_hit.z()) <
// max_dist);
// });
for (const auto& h : hits) {
if(std::hypot(h.x() - ref_hit.x(), h.y() - ref_hit.y(), h.z() - ref_hit.z()) < max_dist){
first_cluster_hits.push_back(h.clone());
} else {
remaining_hits.push_back(h.clone());
}
}
debug() << " first cluster size " << first_cluster_hits.size() << endmsg;
double total_energy = std::accumulate(std::begin(first_cluster_hits), std::end(first_cluster_hits), 0.0,
[](double t,const eic::CalorimeterHit& h1 ) { return (t + h1.energy()); });
debug() << " total_energy = " << total_energy << endmsg;
eic::Cluster cluster0;
//debug() << " cluster0 = " << cluster0 << endmsg;
//eic::Cluster cluster1 = std::accumulate(std::begin(first_cluster_hits), std::end(first_cluster_hits), cluster0,
for (const auto& h : first_cluster_hits) {
cluster0.energy(cluster0.energy() + h.energy());
cluster0.x(cluster0.x() + h.energy() * h.x()/total_energy);
cluster0.y(cluster0.y() + h.energy() * h.y()/total_energy);
cluster0.z(cluster0.z() + h.energy() * h.z()/total_energy);
}
debug() << " cluster0 = " << cluster0 << endmsg;
clusters.push_back(cluster0);
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(SimpleClustering)
} // namespace Jug::Reco
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment