Commit 544e9ea5 authored by Chao Peng's avatar Chao Peng
Browse files

Merge branch 'master' into 'master'

Use local positions for clustering, and convert to global position in reconstruction

See merge request !11
parents 67d2e782 a803b500
......@@ -16,14 +16,14 @@
namespace Jug {
namespace Digi {
/** Crystal Endcaps Calorimeter detector digitization.
*
*
*/
class CrystalEndcapsDigi : public GaudiAlgorithm {
public:
Gaudi::Property<double> m_energyResolution{this, "energyResolution", 0.02}; // 2%sqrt(E)
Rndm::Numbers m_gaussDist;
DataHandle<dd4pod::CalorimeterHitCollection> m_inputHitCollection{ "inputHitCollection", Gaudi::DataHandle::Reader, this};
......@@ -51,11 +51,12 @@ namespace Jug {
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
for(const auto& ahit : *simhits){
eic::RawCalorimeterHit rawhit((long long)ahit.cellID(), (long long)ahit.cellID(),
(long long)(ahit.energyDeposit() + m_gaussDist*sqrt(ahit.energyDeposit()))/Gaudi::Units::MeV * 100.0,
(double)ahit.truth().time);
rawhits->push_back(rawhit);
for (const auto& ahit : *simhits) {
eic::RawCalorimeterHit rawhit(
(long long) ahit.cellID(),
(long long) (ahit.energyDeposit() + m_gaussDist*sqrt(ahit.energyDeposit()))/Gaudi::Units::MeV * 100.0,
(double) ahit.truth().time/Gaudi::Units::ns);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
}
......
......@@ -111,7 +111,7 @@ namespace Jug {
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
for(const auto& ahit : *simhits) {
//std::cout << ahit << "\n";
eic::RawCalorimeterHit rawhit((long long)ahit.cellID(), (long long)ahit.cellID(),
eic::RawCalorimeterHit rawhit((long long)ahit.cellID(),
(long long)ahit.energyDeposit() * 100, 0);
rawhits->push_back(rawhit);
}
......
......@@ -74,7 +74,7 @@ namespace Jug {
double aterm = m_gaussDist_a()*sqrtE;
double bterm = ahit.energyDeposit()*m_gaussDist_b();
// here 1000 is arbitrary scale factor
eic::RawCalorimeterHit rawhit((long long)ahit.cellID(), (long long)ahit.cellID(),
eic::RawCalorimeterHit rawhit((long long)ahit.cellID(),
(long long)(ahit.energyDeposit() +aterm + bterm) * 1000, 0);
rawhits->push_back(rawhit);
}
......
......@@ -42,15 +42,16 @@ public:
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ClusterCollection>
m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer, this};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
DataHandle<eic::CalorimeterHitCollection>
m_splitHitCollection{"splitHitCollection", Gaudi::DataHandle::Writer, this};
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CalorimeterIslandCluster(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputClusterCollection", m_outputClusterCollection, "");
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputClusterCollection", m_outputClusterCollection, "");
declareProperty("splitHitCollection", m_splitHitCollection, "");
}
StatusCode initialize() override
......@@ -58,12 +59,6 @@ public:
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;
}
......@@ -73,26 +68,28 @@ public:
const auto &hits = *m_inputHitCollection.get();
// Create output collections
auto &clusters = *m_outputClusterCollection.createAndPut();
auto &split_hits = *m_splitHitCollection.createAndPut();
// group neighboring hits
std::vector<bool> visits(hits.size(), false);
eic::ClusterCollection groups;
std::vector<std::vector<eic::CalorimeterHit>> groups;
for (size_t i = 0; i < hits.size(); ++i)
{
// already in a group
if (visits[i]) {
continue;
}
groups.emplace_back();
// create a new group, and group all the neighboring hits
dfs_group(groups.create(), i, hits, visits);
dfs_group(groups.back(), i, hits, visits);
}
info() << "we have " << groups.size() << " groups of hits" << endmsg;
// info() << "we have " << groups.size() << " groups of hits" << endmsg;
for (auto &group : groups) {
auto maxima = find_local_maxima(group);
auto split_collection = split_group(group, maxima, clusters);
info() << "hits in a group: " << group.hits_size() << ", "
<< "local maxima: " << maxima.hits_size() << endmsg;
split_group(group, maxima, clusters, split_hits);
// info() << "hits in a group: " << group.hits_size() << ", "
// << "local maxima: " << maxima.hits_size() << endmsg;
}
return StatusCode::SUCCESS;
......@@ -102,23 +99,16 @@ private:
// helper function to group hits
inline bool is_neighbor(const eic::ConstCalorimeterHit &h1, const eic::ConstCalorimeterHit &h2)
{
auto pos1 = h1.position();
auto pos2 = h2.position();
auto dim1 = m_geoSvc->cellIDPositionConverter()->cellDimensions(h1.cellID0());
auto dim2 = m_geoSvc->cellIDPositionConverter()->cellDimensions(h2.cellID0());
// info() << std::abs(pos1.x - pos2.x) << ", " << (dim1[0] + dim2[0])/2. << ", "
// << std::abs(pos1.y - pos2.y) << ", " << (dim1[1] + dim2[1])/2. << endmsg;
return (std::abs(pos1.x - pos2.x) <= (dim1[0] + dim2[0])/2.*m_groupRange) &&
(std::abs(pos1.y - pos2.y) <= (dim1[1] + dim2[1])/2.*m_groupRange);
return (std::abs(h1.local_x() - h2.local_x()) <= (h1.dim_x() + h2.dim_y())/2.*m_groupRange) &&
(std::abs(h1.local_y() - h2.local_y()) <= (h1.dim_y() + h2.dim_y())/2.*m_groupRange);
}
// grouping function with Depth-First Search
void dfs_group(eic::Cluster group, int idx, const eic::CalorimeterHitCollection &hits, std::vector<bool> &visits)
void dfs_group(std::vector<eic::CalorimeterHit> &group, int idx,
const eic::CalorimeterHitCollection &hits, std::vector<bool> &visits)
{
auto hit = hits[idx];
group.addhits(hit);
group.push_back(hit);
visits[idx] = true;
for(size_t i = 0; i < hits.size(); ++i)
{
......@@ -130,10 +120,10 @@ private:
}
// find local maxima that above a certain threshold
eic::Cluster find_local_maxima(const eic::Cluster &group)
eic::Cluster find_local_maxima(const std::vector<eic::CalorimeterHit> &group)
{
eic::Cluster maxima;
for(auto &hit : group.hits())
for(auto &hit : group)
{
// not a qualified center
if(hit.energy() < m_minClusterCenterEdep) {
......@@ -141,7 +131,7 @@ private:
}
bool maximum = true;
for(auto &hit2 : group.hits())
for(auto &hit2 : group)
{
if(hit == hit2)
continue;
......@@ -168,36 +158,34 @@ private:
}
// split a group of hits according to the local maxima
eic::CalorimeterHitCollection split_group(eic::Cluster group, const eic::Cluster &maxima,
eic::ClusterCollection &clusters)
// split_hits is used to persistify the data
void split_group(const std::vector<eic::CalorimeterHit> &group, const eic::Cluster &maxima,
eic::ClusterCollection &clusters, eic::CalorimeterHitCollection &split_hits)
{
// to persistify the split hits
eic::CalorimeterHitCollection scoll;
// special cases
if (maxima.hits_size() == 0) {
return scoll;
return;
} else if (maxima.hits_size() == 1) {
clusters.push_back(group.clone());
return scoll;
auto cl = clusters.create();
for (auto &hit : group) {
cl.addhits(hit);
}
return;
}
// distance reference
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(maxima.hits_begin()->cellID0());
double dist_ref = dim[0];
// split between maxima
std::vector<double> weights(maxima.hits_size());
std::vector<eic::Cluster> splits(maxima.hits_size());
size_t i = 0;
for (auto it = group.hits_begin(); it != group.hits_end(); ++it, ++i) {
auto hpos = it->position();
for (auto it = group.begin(); it != group.end(); ++it, ++i) {
auto hedep = it->energy();
size_t j = 0;
// calculate weights for local maxima
for (auto cit = maxima.hits_begin(); cit != maxima.hits_end(); ++cit, ++j) {
double dist_ref = cit->dim_x();
double energy = cit->energy();
auto pos = cit->position();
double dist = std::sqrt(std::pow(pos.x - hpos.x, 2) + std::pow(pos.y - hpos.y, 2));
double dist = std::sqrt(std::pow(it->local_x() - cit->local_x(), 2)
+ std::pow(it->local_y() - cit->local_y(), 2));
weights[j] = std::exp(-dist/dist_ref)*energy;
}
......@@ -216,10 +204,10 @@ private:
if (weight <= 1e-6) {
continue;
}
eic::CalorimeterHit hit(it->cellID0(), it->cellID1(), hedep*weight,
it->time(), it->position(), it->type());
scoll.push_back(hit);
auto hit = it->clone();
hit.energy(hedep*weight);
hit.type(1);
split_hits.push_back(hit);
splits[k].addhits(hit);
}
}
......@@ -227,7 +215,7 @@ private:
for (auto &cl : splits) {
clusters.push_back(cl);
}
return scoll;
return;
}
};
......
......@@ -34,6 +34,8 @@ public:
Gaudi::Property<double> m_logWeightBase{this, "logWeightBase", 3.6};
DataHandle<eic::ClusterCollection>
m_clusterCollection{"clusterCollection", Gaudi::DataHandle::Reader, this};
// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc)
......@@ -47,42 +49,72 @@ public:
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 &clusters = *m_clusterCollection.get();
auto &clusters = *m_clusterCollection.get();
// reconstruct hit position for the cluster
for (auto &cl : clusters) {
reconstruct(cl);
info() << cl.energy()/GeV << " GeV, (" << cl.position()[0]/mm << ", "
<< cl.position()[1]/mm << ", " << cl.position()[2]/mm << ")" << endmsg;
// info() << cl.energy()/GeV << " GeV, (" << cl.position().x/mm << ", "
// << cl.position().y/mm << ", " << cl.position().z/mm << ")" << endmsg;
}
return StatusCode::SUCCESS;
}
private:
void reconstruct(eic::Cluster cl) {
float totalE = 0.;
void reconstruct(eic::Cluster cl)
{
// no hits
if (cl.hits_size() == 0) {
return;
}
// calculate total energy, find the cell with the maximum energy deposit
float totalE = 0., maxE = 0.;
auto centerID = cl.hits_begin()->cellID();
for (auto &hit : cl.hits()) {
totalE += hit.energy();
auto energy = hit.energy();
totalE += energy;
if (energy > maxE) {
maxE = energy;
centerID = hit.cellID();
}
}
cl.energy(totalE);
// center of gravity with logarithmic weighting
float totalW = 0., x = 0., y = 0., z = 0.;
float tw = 0., x = 0., y = 0., z = 0.;
for (auto &hit : cl.hits()) {
// suppress low energy contributions
float weight = std::max(0., m_logWeightBase + std::log(hit.energy()/totalE));
totalW += weight;
x += hit.position().x * weight;
y += hit.position().y * weight;
z += hit.position().z * weight;
float w = std::max(0., m_logWeightBase + std::log(hit.energy()/totalE));
tw += w;
x += hit.local_x() * w;
y += hit.local_y() * w;
z += hit.local_z() * w;
}
cl.position() = std::array<float, 3>{x/totalW, y/totalW, z/totalW};
// convert local position to global position, use the cell with max edep as a reference
auto volman = m_geoSvc->detector()->volumeManager();
auto alignment = volman.lookupDetector(centerID).nominal();
// depth
// @TODO, assume on the surface
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(centerID);
double depth = -dim[2]/2.;
// info() << depth << " (" << dim[0] << ", " << dim[1] << ", " << dim[2] << ")" << endmsg;
auto gpos = alignment.localToWorld(dd4hep::Position(x/tw, y/tw, z/tw + depth));
cl.position({gpos.x(), gpos.y(), gpos.z()});
}
};
......
......@@ -67,10 +67,20 @@ public:
for (auto &rh : rawhits) {
float energy = rh.amplitude()/100.*MeV;
if (energy >= m_minModuleEdep) {
float time = rh.timeStamp();
auto pos = m_geoSvc->cellIDPositionConverter()->position(rh.cellID0());
float time = rh.timeStamp()*ns;
auto id = rh.cellID();
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
// cell dimension
auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
hits.push_back(eic::CalorimeterHit{
rh.cellID0(), rh.cellID1(), energy, time, {pos.X(), pos.Y(), pos.Z()}, 0
id, energy, time,
{gpos.x(), gpos.y(), gpos.z()},
{pos.x(), pos.y(), pos.z()},
{dim[0], dim[1], dim[2]},
0
});
}
}
......
Markdown is supported
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