Commit 1a24a7e7 authored by Chao Peng's avatar Chao Peng
Browse files

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

parent 4e9e4a73
......@@ -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,13 @@ 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.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;
}
......
......@@ -42,7 +42,7 @@ public:
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ClusterCollection>
m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer, this};
/// Pointer to the geometry service
// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
......@@ -86,13 +86,13 @@ public:
// create a new group, and group all the neighboring hits
dfs_group(groups.create(), 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;
// info() << "hits in a group: " << group.hits_size() << ", "
// << "local maxima: " << maxima.hits_size() << endmsg;
}
return StatusCode::SUCCESS;
......@@ -102,8 +102,9 @@ 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();
// check neighbor hits with local positions
auto pos1 = h1.localPosition();
auto pos2 = h2.localPosition();
auto dim1 = m_geoSvc->cellIDPositionConverter()->cellDimensions(h1.cellID0());
auto dim2 = m_geoSvc->cellIDPositionConverter()->cellDimensions(h2.cellID0());
......@@ -190,13 +191,13 @@ private:
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();
auto hpos = it->localPosition();
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 energy = cit->energy();
auto pos = cit->position();
auto pos = cit->localPosition();
double dist = std::sqrt(std::pow(pos.x - hpos.x, 2) + std::pow(pos.y - hpos.y, 2));
weights[j] = std::exp(-dist/dist_ref)*energy;
}
......@@ -218,7 +219,7 @@ private:
}
eic::CalorimeterHit hit(it->cellID0(), it->cellID1(), hedep*weight,
it->time(), it->position(), it->type());
it->time(), it->localPosition(), it->type());
scoll.push_back(hit);
splits[k].addhits(hit);
}
......
......@@ -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,6 +49,12 @@ 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;
}
......@@ -57,32 +65,51 @@ public:
// 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()->cellID0();
for (auto &hit : cl.hits()) {
totalE += hit.energy();
auto energy = hit.energy();
totalE += energy;
if (energy > maxE) {
maxE = energy;
centerID = hit.cellID0();
}
}
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.localPosition().x * w;
y += hit.localPosition().y * w;
z += hit.localPosition().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();
auto gpos = alignment.localToWorld(dd4hep::Position(x/tw, y/tw, z/tw));
cl.position({gpos.x(), gpos.y(), gpos.z()});
}
};
......
......@@ -67,10 +67,11 @@ 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;
// local positions
auto pos = m_geoSvc->cellIDPositionConverter()->findContext(rh.cellID0())->volumePlacement().position();
hits.push_back(eic::CalorimeterHit{
rh.cellID0(), rh.cellID1(), energy, time, {pos.X(), pos.Y(), pos.Z()}, 0
rh.cellID0(), rh.cellID1(), energy, time, {pos.x(), pos.y(), pos.z()}, 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