Commit ba64dee3 authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

Fix reconstruction

parent 3aac4802
......@@ -85,15 +85,14 @@ juggler:default:
needs:
- version
script:
- CACHEBUST=`date +%s`
- git clone https://eicweb.phy.anl.gov/containers/eic_container.git
- ./eic_container/gitlab-ci/docker_login.sh --ci -n $DOCKER_NTRIES -t $DOCKER_WAIT_TIME
- echo ${CI_COMMIT_BRANCH}
- docker build -t ${CI_REGISTRY_IMAGE}/${BUILD_IMAGE}:${INTERNAL_TAG}
-f eic_container/containers/jug/xl.Dockerfile
--build-arg INTERNAL_TAG=${EIC_DEV_TAG}
--build-arg CACHEBUST=${CACHEBUST}
--build-arg JUGGLER_VERSION=${CI_COMMIT_REF_NAME}
--build-arg JUG_VERSION=juggler-${INTERNAL_TAG}-$(git rev-parse HEAD)
eic_container/containers/jug
- |
if [ "$CI_COMMIT_BRANCH" == "master" ]; then
......
......@@ -11,6 +11,18 @@ find_package(NPDet REQUIRED)
set(PODIO $ENV{PODIO})
set(CMAKE_MODULE_PATH CMAKE_MODULE_PATH PODIO)
find_package(podio 0.11.0 REQUIRED)
## There was a breaking change in how collections are stored in podio between
## version v0.13.0 and v0.13.1. This checks for podio version and sets the
## PODIO_BEFORE_0_13_1 preprocessor variable
if (${podio_VERSION_MAJOR} EQUAL 0 AND ${podio_VERSION_MINOR} LESS_EQUAL 13)
if (${podio_VERSION_MINOR} EQUAL 13 AND ${podio_VERSION_PATCH} EQUAL 0)
add_definitions("-DPODIO_BEFORE_0_13_1")
elseif (${podio_VERSION_MINOR} LESS 13)
add_definitions("-DPODIO_BEFORE_0_13_1")
endif()
endif()
include_directories(${podio_INCLUDE_DIR})
find_package(ROOT COMPONENTS MathCore GenVector Geom REQUIRED)
find_package(DD4hep COMPONENTS DDG4 DDG4IO DDRec REQUIRED)
......
......@@ -2,10 +2,10 @@
// Author: Chao Peng
// Date: 06/14/2021
#include <algorithm>
#include <bitset>
#include "fmt/format.h"
#include "fmt/ranges.h"
#include <algorithm>
#include <bitset>
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
......@@ -31,14 +31,14 @@ using namespace Gaudi::Units;
namespace Jug::Reco {
class CalorimeterHitReco : public GaudiAlgorithm {
public:
class CalorimeterHitReco : public GaudiAlgorithm {
public:
// length unit from dd4hep, should be fixed
Gaudi::Property<double> m_lUnit{this, "lengthUnit", dd4hep::mm};
// digitization settings, must be consistent with digi class
Gaudi::Property<int> m_capADC{this, "capacityADC", 8096};
Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100.*MeV};
Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100. * MeV};
Gaudi::Property<int> m_pedMeanADC{this, "pedestalMean", 400};
Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
......@@ -46,150 +46,164 @@ public:
Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
Gaudi::Property<double> m_minEdep{this, "minimumHitEdep", 0.};
DataHandle<eic::RawCalorimeterHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::CalorimeterHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{
"outputHitCollection", Gaudi::DataHandle::Writer, this};
// geometry service to get ids, ignored if no names provided
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
Gaudi::Property<std::string> m_layerField{this, "layerField", ""};
Gaudi::Property<std::string> m_sectorField{this, "sectorField", ""};
SmartIF<IGeoSvc> m_geoSvc;
dd4hep::BitFieldCoder *id_dec = nullptr;
size_t sector_idx, layer_idx;
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
Gaudi::Property<std::string> m_layerField{this, "layerField", ""};
Gaudi::Property<std::string> m_sectorField{this, "sectorField", ""};
SmartIF<IGeoSvc> m_geoSvc;
dd4hep::BitFieldCoder* id_dec = nullptr;
size_t sector_idx, layer_idx;
// name of detelment or fields to find the local detector (for global->local transform)
// if nothing is provided, the lowest level DetElement (from cellID) will be used
Gaudi::Property<std::string> m_localDetElement{this, "localDetElement", ""};
Gaudi::Property<std::string> m_localDetElement{this, "localDetElement", ""};
Gaudi::Property<std::vector<std::string>> u_localDetFields{this, "localDetFields", {}};
dd4hep::Alignment local;
size_t local_mask = ~0;
dd4hep::Alignment local;
size_t local_mask = ~0;
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CalorimeterHitReco(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
m_geoSvc = service(m_geoSvcName);
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;
}
// do not get the layer/sector ID if no readout class provided
if (m_readout.value().empty()) {
return StatusCode::SUCCESS;
}
m_geoSvc = service(m_geoSvcName);
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;
auto id_spec = m_geoSvc->detector()->readout(m_readout).idSpec();
try {
id_dec = id_spec.decoder();
if (m_sectorField.value().size()) {
sector_idx = id_dec->index(m_sectorField);
}
// do not get the layer/sector ID if no readout class provided
if (m_readout.value().empty()) {
return StatusCode::SUCCESS;
if (m_layerField.value().size()) {
layer_idx = id_dec->index(m_layerField);
}
} catch (...) {
error() << "Failed to load ID decoder for " << m_readout << endmsg;
return StatusCode::FAILURE;
}
auto id_spec = m_geoSvc->detector()->readout(m_readout).idSpec();
// local detector name has higher priority
if (m_localDetElement.value().size()) {
try {
id_dec = id_spec.decoder();
if (m_sectorField.value().size()) { sector_idx = id_dec->index(m_sectorField); }
if (m_layerField.value().size()) { layer_idx = id_dec->index(m_layerField); }
local = m_geoSvc->detector()->detector(m_localDetElement.value()).nominal();
info() << "Local coordinate system from DetElement " << m_localDetElement.value()
<< endmsg;
} catch (...) {
error() << "Failed to load ID decoder for " << m_readout << endmsg;
return StatusCode::FAILURE;
error() << "Failed to locate local coordinate system from DetElement "
<< m_localDetElement.value() << endmsg;
return StatusCode::FAILURE;
}
// local detector name has higher priority
if (m_localDetElement.value().size()) {
try {
local = m_geoSvc->detector()->detector(m_localDetElement.value()).nominal();
info() << "Local coordinate system from DetElement " << m_localDetElement.value() << endmsg;
} catch (...) {
error() << "Failed to locate local coordinate system from DetElement " << m_localDetElement.value()
<< endmsg;
return StatusCode::FAILURE;
}
// or get from fields
} else {
std::vector<std::pair<std::string, int>> fields;
for (auto &f : u_localDetFields.value()) {
fields.push_back({f, 0});
}
local_mask = id_spec.get_mask(fields);
// use all fields if nothing provided
if (fields.empty()) { local_mask = ~0; }
info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask, fmt::join(fields, ", "))
<< endmsg;
} else {
std::vector<std::pair<std::string, int>> fields;
for (auto& f : u_localDetFields.value()) {
fields.push_back({f, 0});
}
local_mask = id_spec.get_mask(fields);
// use all fields if nothing provided
if (fields.empty()) {
local_mask = ~0;
}
info() << fmt::format("Local DetElement mask {:#064b} from fields [{}]", local_mask,
fmt::join(fields, ", "))
<< endmsg;
}
return StatusCode::SUCCESS;
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collections
const auto& rawhits = *m_inputHitCollection.get();
// create output collections
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
for (auto& rh : rawhits) {
// did not pass the zero-suppression threshold
if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC*m_pedSigmaADC) {
continue;
}
// convert ADC -> energy
float energy = (rh.amplitude() - m_pedMeanADC) / (float) m_capADC * m_dyRangeADC;
if (energy < m_minEdep) {
continue;
}
float time = rh.timeStamp(); // ns
auto id = rh.cellID();
int lid = ((id_dec != nullptr) & m_layerField.value().size())
? static_cast<int>(id_dec->get(id, layer_idx))
: -1;
int sid = ((id_dec != nullptr) & m_sectorField.value().size())
? static_cast<int>(id_dec->get(id, sector_idx))
: -1;
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
if (m_localDetElement.value().empty()) {
auto volman = m_geoSvc->detector()->volumeManager();
local = volman.lookupDetElement(id & local_mask).nominal();
}
auto pos = local.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
// auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
// cell dimension
auto cdim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
double dim[3] = {0., 0., 0.};
for (size_t i = 0; i < cdim.size() && i < 3; ++i) {
dim[i] = cdim[i] / m_lUnit;
}
// debug() << std::bitset<64>(id) << "\n"
// << m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().volIDs().str()
// << endmsg;
hits.push_back(eic::CalorimeterHit{
id, -1, lid, sid, // cell id, cluster id, layer id, sector id
energy, time, // energy, time
{gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit},
{pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit},
{dim[0], dim[1], dim[2]},
0 // @TODO: hit type
});
// input collections
const auto& rawhits = *m_inputHitCollection.get();
// create output collections
auto& hits = *m_outputHitCollection.createAndPut();
// energy time reconstruction
for (const auto& rh : rawhits) {
// did not pass the zero-suppression threshold
if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC * m_pedSigmaADC) {
continue;
}
return StatusCode::SUCCESS;
// convert ADC -> energy
float energy = (rh.amplitude() - m_pedMeanADC) / (float)m_capADC * m_dyRangeADC;
if (energy < m_minEdep) {
continue;
}
float time = rh.timeStamp(); // ns
auto id = rh.cellID();
int lid = ((id_dec != nullptr) & m_layerField.value().size())
? static_cast<int>(id_dec->get(id, layer_idx))
: -1;
int sid = ((id_dec != nullptr) & m_sectorField.value().size())
? static_cast<int>(id_dec->get(id, sector_idx))
: -1;
// global positions
auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
// local positions
if (m_localDetElement.value().empty()) {
auto volman = m_geoSvc->detector()->volumeManager();
local = volman.lookupDetElement(id & local_mask).nominal();
}
auto pos = local.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
// auto pos =
// m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position(); cell
// dimension
auto cdim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
double dim[3] = {0., 0., 0.};
for (size_t i = 0; i < cdim.size() && i < 3; ++i) {
dim[i] = cdim[i] / m_lUnit;
}
// debug() << std::bitset<64>(id) << "\n"
// <<
// m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().volIDs().str()
// << endmsg;
hits.push_back({
id,
-1,
lid,
sid, // cell id, cluster id, layer id, sector id
energy,
time, // energy, time
{gpos.x() / m_lUnit, gpos.y() / m_lUnit, gpos.z() / m_lUnit},
{pos.x() / m_lUnit, pos.y() / m_lUnit, pos.z() / m_lUnit},
{dim[0], dim[1], dim[2]},
0 // @TODO: hit type
});
}
return StatusCode::SUCCESS;
}
}; // class CalorimeterHitReco
}; // class CalorimeterHitReco
DECLARE_COMPONENT(CalorimeterHitReco)
DECLARE_COMPONENT(CalorimeterHitReco)
} // namespace Jug::Reco
......@@ -4,25 +4,25 @@
*
* Author: Chao Peng (ANL), 03/31/2021
*/
#include <bitset>
#include <algorithm>
#include <bitset>
#include <unordered_map>
#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 "GaudiAlg/Transformer.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/ToolHandle.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DDRec/SurfaceManager.h"
#include "DDSegmentation/BitFieldCoder.h"
#include "fmt/ranges.h"
#include "fmt/format.h"
#include "fmt/ranges.h"
// FCCSW
#include "JugBase/DataHandle.h"
......@@ -35,124 +35,124 @@ using namespace Gaudi::Units;
namespace Jug::Reco {
class CalorimeterHitsMerger : public GaudiAlgorithm {
public:
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
class CalorimeterHitsMerger : public GaudiAlgorithm {
public:
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
// field names to generate id mask, the hits will be grouped by masking the field
Gaudi::Property<std::vector<std::string>> u_fields{this, "fields", {"layer"}};
Gaudi::Property<std::vector<std::string>> u_fields{this, "fields", {"layer"}};
// reference field numbers to locate position for each merged hits group
Gaudi::Property<std::vector<int>> u_refs{this, "fieldRefNumbers", {}};
DataHandle<eic::CalorimeterHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::CalorimeterHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
Gaudi::Property<std::vector<int>> u_refs{this, "fieldRefNumbers", {}};
DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection",
Gaudi::DataHandle::Reader, this};
DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{
"outputHitCollection", Gaudi::DataHandle::Writer, this};
SmartIF<IGeoSvc> m_geoSvc;
uint64_t id_mask, ref_mask;
uint64_t id_mask, ref_mask;
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CalorimeterHitsMerger(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
m_geoSvc = service(m_geoSvcName);
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;
}
if (m_readout.value().empty()) {
error() << "readoutClass is not provided, it is needed to know the fields in readout ids" << endmsg;
return StatusCode::FAILURE;
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
m_geoSvc = service(m_geoSvcName);
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;
}
if (m_readout.value().empty()) {
error() << "readoutClass is not provided, it is needed to know the fields in readout ids"
<< endmsg;
return StatusCode::FAILURE;
}
try {
auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec();
id_mask = 0;
std::vector<std::pair<std::string, int>> ref_fields;
for (size_t i = 0; i < u_fields.size(); ++i) {
id_mask |= id_desc.field(u_fields[i])->mask();
// use the provided id number to find ref cell, or use 0
int ref = i < u_refs.size() ? u_refs[i] : 0;
ref_fields.push_back({u_fields[i], ref});
}
try {
auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec();
id_mask = 0;
std::vector<std::pair<std::string, int>> ref_fields;
for (size_t i = 0; i < u_fields.size(); ++i) {
id_mask |= id_desc.field(u_fields[i])->mask();
// use the provided id number to find ref cell, or use 0
int ref = i < u_refs.size() ? u_refs[i] : 0;
ref_fields.push_back({u_fields[i], ref});
}
ref_mask = id_desc.encode(ref_fields);
// debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg;
} catch (...) {
error() << "Failed to load ID decoder for " << m_readout << endmsg;
return StatusCode::FAILURE;
}
id_mask = ~id_mask;
info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask)
<< endmsg;
return StatusCode::SUCCESS;
ref_mask = id_desc.encode(ref_fields);
// debug() << fmt::format("Referece id mask for the fields {:#064b}", ref_mask) << endmsg;
} catch (...) {
error() << "Failed to load ID decoder for " << m_readout << endmsg;
return StatusCode::FAILURE;
}
id_mask = ~id_mask;
info() << fmt::format("ID mask in {:s}: {:#064b}", m_readout, id_mask) << endmsg;
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collections
const auto &hits = *m_inputHitCollection.get();
// Create output collections
auto &mhits = *m_outputHitCollection.createAndPut();
// dd4hep decoders
auto poscon = m_geoSvc->cellIDPositionConverter();
auto volman = m_geoSvc->detector()->volumeManager();
// sum energies that has the same id
std::unordered_map<long long, size_t> merge_map;
for (auto &h : hits) {
int64_t id = h.cellID() & id_mask;
// use the reference field position
int64_t ref_id = id | ref_mask;
// debug() << h.cellID() << " - " << std::bitset<64>(h.cellID()) << endmsg;
auto it = merge_map.find(id);
debug() << fmt::format("{:#064b} - {:#064b}, ref: {:#064b}",
h.cellID(), id, ref_id) << endmsg;
if (it == merge_map.end()) {
merge_map[id] = mhits.size();
auto ahit = h.clone();
ahit.cellID(ref_id);
// global positions
auto gpos = poscon->position(ref_id);
// local positions
auto alignment = volman.lookupDetElement(ref_id).nominal();
// debug() << volman.lookupDetElement(ref_id).path() << ", " << volman.lookupDetector(ref_id).path() << endmsg;
auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
ahit.position({gpos.x()/dd4hep::mm, gpos.y()/dd4hep::mm, gpos.z()/dd4hep::mm});
ahit.local({pos.x()/dd4hep::mm, pos.y()/dd4hep::mm, pos.z()/dd4hep::mm});
mhits.push_back(ahit);
} else {
mhits[it->second].energy(mhits[it->second].energy() + h.energy());
}
// input collections
const auto& hits = *m_inputHitCollection.get();
// Create output collections
auto& mhits = *m_outputHitCollection.createAndPut();
// dd4hep decoders
auto poscon = m_geoSvc->cellIDPositionConverter();
auto volman = m_geoSvc->detector()->volumeManager();
// sum energies that has the same id
std::unordered_map<long long, size_t> merge_map;
for (const auto& h : hits) {
int64_t id = h.cellID() & id_mask;
// use the reference field position
int64_t ref_id = id | ref_mask;
// debug() << h.cellID() << " - " << std::bitset<64>(h.cellID()) << endmsg;
auto it = merge_map.find(id);
debug() << fmt::format("{:#064b} - {:#064b}, ref: {:#064b}", h.cellID(), id, ref_id)
<< endmsg;
if (it == merge_map.end()) {
merge_map[id] = mhits.size();
// global positions
auto gpos = poscon->position(ref_id);
// local positions
auto alignment = volman.lookupDetElement(ref_id).nominal();
// debug() << volman.lookupDetElement(ref_id).path() << ", " <<
// volman.lookupDetector(ref_id).path() << endmsg;
auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
mhits.push_back({ref_id,
h.clusterID(),
h.layerID(),
h.sectorID(),
h.energy(),
h.time(),
{gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm},
{pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm},
h.dimension(),
h.type()});
} else {
mhits[it->second].energy(mhits[it->second].energy() + h.energy());
}
}