Commit 19b252d9 authored by Jihee Kim's avatar Jihee Kim Committed by Chao Peng
Browse files

Resolve "Merger Algorithm"

parent 8ef2bef5
......@@ -47,6 +47,7 @@ gaudi_add_module(JugRecoPlugins
src/components/ImagingPixelMerger.cpp
src/components/ImagingTopoCluster.cpp
src/components/ImagingClusterReco.cpp
src/components/CalorimeterHitsEtaPhiProjector.cpp
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec )
target_compile_options(JugRecoPlugins PRIVATE -Wno-suggest-override)
......
/*
* A hits converter to prepare dataset for machine learning
* It converts hits with (x, y, z, E) to (E, eta, phi) layer by layer
* With a defined grid size and ranges, it merge the hits within one grid and drop-off hits
* out-of-range The capacity of each layer is fixed (padding with zeros), and the hits with least
* energies that exceed the capacity will be discarded.
*
* Author: Chao Peng (ANL), Jihee Kim, 07/16/2021
*
*/
#include <algorithm>
#include <bitset>
#include <unordered_map>
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/ToolHandle.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/Surface.h"
#include "DDRec/SurfaceManager.h"
// FCCSW
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugBase/Utilities/Utils.hpp"
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
using namespace Gaudi::Units;
typedef ROOT::Math::XYZPoint Point3D;
struct pair_hash {
template <class T1, class T2>
std::size_t operator()(const std::pair<T1, T2>& pair) const
{
return std::hash<T1>()(pair.first) ^ std::hash<T2>()(pair.second);
}
};
namespace Jug::Reco {
class CalorimeterHitsEtaPhiProjector : public GaudiAlgorithm {
public:
Gaudi::Property<std::vector<double>> u_gridSizes{this, "gridSizes", {0.001, 0.001*rad}};
DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{
"outputHitCollection", Gaudi::DataHandle::Writer, this};
double gridSizes[2];
CalorimeterHitsEtaPhiProjector(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
if (u_gridSizes.size() != 2) {
error() << "Expected 2 values for gridSizes, received " << u_gridSizes.size() << endmsg;
return StatusCode::FAILURE;
}
gridSizes[0] = u_gridSizes.value()[0];
gridSizes[1] = u_gridSizes.value()[1] / rad;
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// Create output collections
auto& mhits = *m_outputHitCollection.createAndPut();
// container
std::unordered_map<std::pair<int64_t, int64_t>, std::vector<eic::ConstCalorimeterHit>, pair_hash> merged_hits;
for (const auto h : *m_inputHitCollection.get()) {
Point3D p(h.x(), h.y(), h.z());
auto bins = std::make_pair(static_cast<int64_t>(pos2bin(p.eta(), gridSizes[0], 0.)),
static_cast<int64_t>(pos2bin(p.phi(), gridSizes[1], 0.)));
merged_hits[bins].push_back(h);
}
for (const auto &[bins, hits] : merged_hits) {
const auto ref = hits.front();
eic::CalorimeterHit hit;
hit.cellID(ref.cellID());
hit.sectorID(ref.sectorID());
hit.layerID(ref.layerID());
// TODO, we can do timing cut to reject noises
hit.time(ref.time());
double r = std::hypot(ref.x(), ref.y(), ref.z());
double eta = bin2pos(bins.first, gridSizes[0], 0.);
double phi = bin2pos(bins.second, gridSizes[1], 1.);
double theta = std::atan(std::exp(-eta))*2.;
ROOT::Math::Polar3DVector p(r, theta, phi);
// TODO, we can do global to local conversion, which gives a better estimate
hit.local(ref.local());
hit.position({p.x(), p.y(), p.z()});
hit.dimension({gridSizes[0], gridSizes[1], 0.});
// merge energy
hit.energy(0.);
for (const auto &h : hits) {
hit.energy(hit.energy() + h.energy());
}
mhits.push_back(hit);
}
return StatusCode::SUCCESS;
}
static int64_t pos2bin(double val, double cell, double offset = 0.) {
return int64_t(std::floor((val + 0.5 * cell - offset) / cell));
}
static double bin2pos(int64_t bin, double cell, double offset = 0.) {
return bin * cell + offset;
}
}; // class CalorimeterHitsEtaPhiProjector
DECLARE_COMPONENT(CalorimeterHitsEtaPhiProjector)
} // namespace Jug::Reco
......@@ -325,6 +325,7 @@ namespace Jug::Reco {
}
// split between maxima
// TODO, here we can implement iterations with profile, or even ML for better splits
std::vector<double> weights(maxima.size());
std::vector<eic::Cluster> splits(maxima.size());
size_t n_clus = clusterID + 1;
......
......@@ -190,9 +190,14 @@ namespace Jug::Reco {
x += hit.x() * w;
y += hit.y() * w;
z += hit.z() * w;
debug() << hit.cellID() << "(" << hit.local_x() << ", " << hit.local_y() << ", "
/*
debug() << hit.cellID() << ": (" << hit.local_x() << ", " << hit.local_y() << ", "
<< hit.local_z() << "), "
<< "(" << hit.x() << ", " << hit.y() << ", " << hit.z() << "), " << endmsg;
*/
}
if (tw == 0.) {
warning() << "zero total weights encountered, you may want to adjust your weighting parameter." << endmsg;
}
res.position({x / tw, y / tw, z / tw});
// convert global position to local position, use the cell with max edep as a reference
......
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