CalorimeterHitDigi.cpp 9.55 KiB
// 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 <algorithm>
#include <cmath>
#include <unordered_map>
#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 "JugBase/IGeoSvc.h"
#include "JugBase/DataHandle.h"
#include "JugBase/UniqueID.h"
#include "fmt/format.h"
#include "fmt/ranges.h"
// Event Model related classes
#include "dd4pod/CalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitData.h"
using namespace Gaudi::Units;
namespace Jug::Digi {
/** Generic calorimeter hit digitiziation.
*
* \ingroup digi
* \ingroup calorimetry
*/
class CalorimeterHitDigi : public GaudiAlgorithm, AlgorithmIDMixin<> {
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};
// digitization settings
Gaudi::Property<int> m_capADC{this, "capacityADC", 8096};
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};
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, stepTDC, tRes, eRes[3] = {0., 0., 0.};
Rndm::Numbers m_normDist;
SmartIF<IGeoSvc> m_geoSvc;
uint64_t id_mask, ref_mask;
DataHandle<dd4pod::CalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawCalorimeterHitCollection> m_outputHitCollection{
"outputHitCollection", Gaudi::DataHandle::Writer, this};
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
CalorimeterHitDigi(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
, AlgorithmIDMixin{name, info()}
{
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().size()) {
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->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;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
if (u_fields.value().size()) {
signal_sum_digi();
} else {
single_hits_digi();
}
return StatusCode::SUCCESS;
}
private:
void single_hits_digi() {
// input collections
const auto simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
int nhits = 0;
for (const auto& ahit : *simhits) {
// Note: juggler internal unit of energy is GeV
// --> optional energy correction to allow for modifying the response for
// effective impelentations, e.g. a homogenous implementation of a fiber
// calorimeter.
const double eDep = ahit.energyDeposit() * m_corrMeanScale;
// apply additional calorimeter noise to corrected energy deposit
const double eResRel = (eDep > 1e-6)
? 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 * (1. + eResRel) / dyRangeADC * m_capADC);
const long long tdc = std::llround((ahit.truth().time + m_normDist() * tRes) * stepTDC);
eic::RawCalorimeterHit rawhit(
{nhits++, algorithmID()},
(long long)ahit.cellID(),
(adc > m_capADC.value() ? m_capADC.value() : adc),
tdc
);
rawhits->push_back(rawhit);
}
}
void signal_sum_digi() {
const auto 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<dd4pod::ConstCalorimeterHit>> merge_map;
for (const auto &ahit : *simhits) {
int64_t hid = (ahit.cellID() & 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
int nhits = 0;
for (auto &[id, hits] : merge_map) {
double edep = hits[0].energyDeposit();
double time = hits[0].truth().time;
double max_edep = hits[0].energyDeposit();
// sum energy, take time from the most energetic hit
for (size_t i = 1; i < hits.size(); ++i) {
edep += hits[i].energyDeposit();
if (hits[i].energyDeposit() > max_edep) {
max_edep = hits[i].energyDeposit();
time = hits[i].truth().time;
}
}
double eResRel = 0.;
// safety check
if (edep > 1e-6) {
eResRel = m_normDist() * eRes[0] / std::sqrt(edep) +
m_normDist() * eRes[1] +
m_normDist() * eRes[2] / edep;
}
double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
long long adc = std::llround(ped + edep * (1. + eResRel) / dyRangeADC * m_capADC);
long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC);
eic::RawCalorimeterHit rawhit(
{nhits++, algorithmID()},
id,
(adc > m_capADC.value() ? m_capADC.value() : adc),
tdc
);
rawhits->push_back(rawhit);
}
}
};
DECLARE_COMPONENT(CalorimeterHitDigi)
} // namespace Jug::Digi