Skip to content

Commit

Permalink
feat: use calorimeter algorithms from EICrecon
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdconinc
wdconinc committed Jun 3, 2024
1 parent cbedab7 commit c4cd4f5
Show file tree
Hide file tree
Showing 12 changed files with 40 additions and 1,982 deletions.
2 changes: 2 additions & 0 deletions JugDigi/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@ gaudi_add_module(JugDigiPlugins
LINK
Gaudi::GaudiAlgLib Gaudi::GaudiKernel
JugBase
JugAlgo
EICrecon::algorithms_calorimetry_library
ROOT::Core ROOT::RIO ROOT::Tree
EDM4HEP::edm4hep
EDM4EIC::edm4eic
Expand Down
264 changes: 5 additions & 259 deletions JugDigi/src/components/CalorimeterHitDigi.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,263 +1,9 @@
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2022 Chao Peng, Wouter Deconinck, Sylvester Joosten
// Copyright (C) 2024 Wouter Deconinck

// A general digitization for CalorimeterHit from simulation
// 1. Smear energy deposit with a/sqrt(E/GeV) + b + c/E or a/sqrt(E/GeV) (relative value)
// 2. Digitize the energy with dynamic ADC range and add pedestal (mean +- sigma)
// 3. Time conversion with smearing resolution (absolute value)
// 4. Signal is summed if the SumFields are provided
//
// Author: Chao Peng
// Date: 06/02/2021
#include <JugAlgo/Algorithm.h>
#include <EICrecon/algorithms/calorimetry/CalorimeterHitDigi.h>

#include <algorithm>
#include <cmath>
#include <unordered_map>
// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
JUGALGO_DEFINE_ALGORITHM(CalorimeterHitDigi, eicrecon::CalorimeterHitDigi, Jug::Digi)

#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "Gaudi/Property.h"
#include "GaudiKernel/RndmGenerators.h"

#include "DDRec/CellIDPositionConverter.h"
#include "DDSegmentation/BitFieldCoder.h"

#include <k4Interface/IGeoSvc.h>
#include <k4FWCore/DataHandle.h>

#include "fmt/format.h"
#include "fmt/ranges.h"

// Event Model related classes
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "edm4eic/RawCalorimeterHitCollection.h"
#include "edm4eic/RawCalorimeterHitData.h"


using namespace Gaudi::Units;

namespace Jug::Digi {

/** Generic calorimeter hit digitiziation.
*
* \ingroup digi
* \ingroup calorimetry
*/
class CalorimeterHitDigi : public GaudiAlgorithm {
public:
// additional smearing resolutions
Gaudi::Property<std::vector<double>> u_eRes{this, "energyResolutions", {}}; // a/sqrt(E/GeV) + b + c/(E/GeV)
Gaudi::Property<double> m_tRes{this, "timeResolution", 0.0 * ns};
// single hit energy deposition threshold
Gaudi::Property<double> m_threshold{this, "threshold", 1. * keV};

// digitization settings
Gaudi::Property<unsigned int> m_capADC{this, "capacityADC", 8096};
Gaudi::Property<double> m_dyRangeADC{this, "dynamicRangeADC", 100 * MeV};
Gaudi::Property<unsigned int> m_pedMeanADC{this, "pedestalMean", 400};
Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
Gaudi::Property<double> m_resolutionTDC{this, "resolutionTDC", 10 * ps};

Gaudi::Property<double> m_corrMeanScale{this, "scaleResponse", 1.0};
// These are better variable names for the "energyResolutions" array which is a bit
// magic @FIXME
//Gaudi::Property<double> m_corrSigmaCoeffE{this, "responseCorrectionSigmaCoeffE", 0.0};
//Gaudi::Property<double> m_corrSigmaCoeffSqrtE{this, "responseCorrectionSigmaCoeffSqrtE", 0.0};

// signal sums
// @TODO: implement signal sums with timing
// field names to generate id mask, the hits will be grouped by masking the field
Gaudi::Property<std::vector<std::string>> u_fields{this, "signalSumFields", {}};
// ref field ids are used for the merged hits, 0 is used if nothing provided
Gaudi::Property<std::vector<int>> u_refs{this, "fieldRefNumbers", {}};
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};

// unitless counterparts of inputs
double dyRangeADC{0}, stepTDC{0}, tRes{0}, eRes[3] = {0., 0., 0.};
Rndm::Numbers m_normDist;
SmartIF<IGeoSvc> m_geoSvc;
uint64_t id_mask{0}, ref_mask{0};

DataHandle<edm4hep::SimCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4eic::RawCalorimeterHitCollection> m_outputHitCollection{
"outputHitCollection", Gaudi::DataHandle::Writer, this};

// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CalorimeterHitDigi(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;
}
// random number generator from service
auto randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
auto sc = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0));
if (!sc.isSuccess()) {
return StatusCode::FAILURE;
}
// set energy resolution numbers
for (size_t i = 0; i < u_eRes.size() && i < 3; ++i) {
eRes[i] = u_eRes[i];
}

// using juggler internal units (GeV, mm, radian, ns)
dyRangeADC = m_dyRangeADC.value() / GeV;
tRes = m_tRes.value() / ns;
stepTDC = ns / m_resolutionTDC.value();

// need signal sum
if (!u_fields.value().empty()) {
m_geoSvc = service(m_geoSvcName);
// sanity checks
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;
}

// get decoders
try {
auto id_desc = m_geoSvc->getDetector()->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.emplace_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.value(), id_mask) << endmsg;
return StatusCode::SUCCESS;
}

return StatusCode::SUCCESS;
}

StatusCode execute() override
{
if (!u_fields.value().empty()) {
signal_sum_digi();
} else {
single_hits_digi();
}
return StatusCode::SUCCESS;
}

private:
void single_hits_digi() {
// input collections
const auto* const simhits = m_inputHitCollection.get();
// Create output collections
auto* rawhits = m_outputHitCollection.createAndPut();
for (const auto& ahit : *simhits) {
// Note: juggler internal unit of energy is GeV
const double eDep = ahit.getEnergy();

// apply additional calorimeter noise to corrected energy deposit
const double eResRel = (eDep > m_threshold)
? m_normDist() * std::sqrt(
std::pow(eRes[0] / std::sqrt(eDep), 2) +
std::pow(eRes[1], 2) +
std::pow(eRes[2] / (eDep), 2)
)
: 0;

const double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / dyRangeADC * m_capADC);

double time = std::numeric_limits<double>::max();
for (const auto& c : ahit.getContributions()) {
if (c.getTime() <= time) {
time = c.getTime();
}
}
const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC);

rawhits->create(
ahit.getCellID(),
(adc > m_capADC.value() ? m_capADC.value() : adc),
tdc
);
}
}

void signal_sum_digi() {
const auto* const simhits = m_inputHitCollection.get();
auto* rawhits = m_outputHitCollection.createAndPut();

// find the hits that belong to the same group (for merging)
std::unordered_map<long long, std::vector<edm4hep::SimCalorimeterHit>> merge_map;
for (const auto &ahit : *simhits) {
int64_t hid = (ahit.getCellID() & id_mask) | ref_mask;
auto it = merge_map.find(hid);

if (it == merge_map.end()) {
merge_map[hid] = {ahit};
} else {
it->second.push_back(ahit);
}
}

// signal sum
for (auto &[id, hits] : merge_map) {
double edep = hits[0].getEnergy();
double time = hits[0].getContributions(0).getTime();
double max_edep = hits[0].getEnergy();
// sum energy, take time from the most energetic hit
for (size_t i = 1; i < hits.size(); ++i) {
edep += hits[i].getEnergy();
if (hits[i].getEnergy() > max_edep) {
max_edep = hits[i].getEnergy();
for (const auto& c : hits[i].getContributions()) {
if (c.getTime() <= time) {
time = c.getTime();
}
}
}
}

// safety check
const double eResRel = (edep > m_threshold)
? m_normDist() * eRes[0] / std::sqrt(edep) +
m_normDist() * eRes[1] +
m_normDist() * eRes[2] / edep
: 0;

double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
unsigned long long adc = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC);
unsigned long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC);

rawhits->create(
id,
(adc > m_capADC.value() ? m_capADC.value() : adc),
tdc
);
}
}
};
// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
DECLARE_COMPONENT(CalorimeterHitDigi)

} // namespace Jug::Digi
1 change: 1 addition & 0 deletions JugFast/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ gaudi_add_module(JugFastPlugins
Gaudi::GaudiAlgLib Gaudi::GaudiKernel
JugBase
JugAlgo
EICrecon::algorithms_calorimetry_library
algorithms::core algorithms::truth
ROOT::Core ROOT::RIO ROOT::Tree
EDM4HEP::edm4hep
Expand Down
88 changes: 4 additions & 84 deletions JugFast/src/components/TruthClustering.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,88 +1,8 @@
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2022 Sylvester Joosten, Whitney Armstrong, Wouter Deconinck
// Copyright (C) 2024 Wouter Deconinck

#include <algorithm>

#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"

#include <k4FWCore/DataHandle.h>
#include <k4Interface/IGeoSvc.h>

// Event Model related classes
#include "edm4hep/MCParticle.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "edm4eic/CalorimeterHitCollection.h"
#include "edm4eic/ClusterCollection.h"
#include "edm4eic/ProtoClusterCollection.h"
#include "edm4eic/RawCalorimeterHitCollection.h"

using namespace Gaudi::Units;

namespace Jug::Fast {

/** Truth clustering algorithm.
*
* \ingroup reco
*/
class TruthClustering : public GaudiAlgorithm {
private:
DataHandle<edm4eic::CalorimeterHitCollection> m_inputHits{"inputHits", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::SimCalorimeterHitCollection> m_mcHits{"mcHits", Gaudi::DataHandle::Reader, this};
DataHandle<edm4eic::ProtoClusterCollection> m_outputProtoClusters{"outputProtoClusters", Gaudi::DataHandle::Writer, this};

public:
TruthClustering(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
declareProperty("inputHits", m_inputHits, "Input calorimeter reco hits");
declareProperty("mcHits", m_mcHits, "Input truth hits");
declareProperty("outputProtoClusters", m_outputProtoClusters, "Output proto clusters");
}

StatusCode initialize() override {
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}

StatusCode execute() override {
// input collections
const auto& hits = *m_inputHits.get();
const auto& mc = *m_mcHits.get();
// Create output collections
auto& proto = *m_outputProtoClusters.createAndPut();

// Map mc track ID to protoCluster index
std::map<int32_t, int32_t> protoIndex;

// Loop over al calorimeter hits and sort per mcparticle
for (const auto& hit : hits) {
const auto& mcHit = mc[hit.getObjectID().index];
const auto& trackID = mcHit.getContributions(0).getParticle().getObjectID().index;
// Create a new protocluster if we don't have one for this trackID
if (protoIndex.count(trackID) == 0) {
auto pcl = proto.create();
protoIndex[trackID] = proto.size() - 1;
}
// Add hit to the appropriate protocluster
proto[protoIndex[trackID]].addToHits(hit);
proto[protoIndex[trackID]].addToWeights(1);
}
return StatusCode::SUCCESS;
}
};
#include <JugAlgo/Algorithm.h>
#include <EICrecon/algorithms/calorimetry/CalorimeterTruthClustering.h>

// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
DECLARE_COMPONENT(TruthClustering)

} // namespace Jug::Fast
JUGALGO_DEFINE_ALGORITHM(CalorimeterTruthClustering, eicrecon::CalorimeterTruthClustering, Jug::Fast)
4 changes: 2 additions & 2 deletions JugReco/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,8 +11,8 @@ gaudi_add_module(JugRecoPlugins
${JugRecoPlugins_sources}
LINK
Gaudi::GaudiAlgLib Gaudi::GaudiKernel
JugBase JugAlgo
algorithms::core algorithms::calorimetry
JugBase
JugAlgo
EICrecon::algorithms_calorimetry_library
ROOT::Core ROOT::RIO ROOT::Tree
EDM4HEP::edm4hep
Expand Down
Loading

0 comments on commit c4cd4f5

Please sign in to comment.