Skip to content

Commit

Permalink
feat: Calorimeter hit digi simplification (#794)
Browse files Browse the repository at this point in the history 
### Briefly, what does this PR introduce?
This applies the #666 treatment to the CalorimeterHitDigi algorithms:
- put all configuration in a single configuration struct for use with
WithPodConfig mixin,
- update the unit tests that are affected,
- create a templated factory-factory to use in the plugins,
- replace all old explicit RawCalorimeterHit factories with the new
factory-factory.

### What kind of change does this PR introduce?
- [ ] Bug fix (issue #__)
- [x] New feature (issue #__)
- [ ] Documentation update
- [ ] Other: __

### Please check if this PR fulfills the following:
- [ ] Tests for the changes have been added
- [ ] Documentation has been added / updated
- [x] Changes have been communicated to collaborators @veprbl
@nathanwbrei

### Does this PR introduce breaking changes? What changes might users
need to make to their code?
No.

### Does this PR change default behavior?
No.

---------

Co-authored-by: Dmitry Kalinkin <[email protected]>
  • Loading branch information
@wdconinc @veprbl
wdconinc and veprbl authored Jul 27, 2023
1 parent 52a912b commit 20c9d80
Show file tree
Hide file tree
Showing 32 changed files with 438 additions and 1,563 deletions.
3 changes: 0 additions & 3 deletions src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,9 +56,6 @@ namespace eicrecon {
const dd4hep::Detector* m_detector;
std::shared_ptr<spdlog::logger> m_log;

std::string m_input_simhit_tag;
std::string m_input_protoclust_tag;

std::function<double(double, double, double, int)> weightFunc;

private:
Expand Down
121 changes: 46 additions & 75 deletions src/algorithms/calorimetry/CalorimeterHitDigi.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,10 +18,8 @@
#include <fmt/format.h>
using namespace dd4hep;

//
// This algorithm converted from:
//
// https://eicweb.phy.anl.gov/EIC/juggler/-/blob/master/JugDigi/src/components/CalorimeterHitDigi.cpp
namespace eicrecon {

//
// TODO:
// - Array type configuration parameters are not yet supported in JANA (needs to be added)
Expand All @@ -31,13 +29,9 @@ using namespace dd4hep;
// values here. This needs to be confirmed.



//------------------------
// AlgorithmInit
//------------------------
void CalorimeterHitDigi::AlgorithmInit(std::shared_ptr<spdlog::logger>& logger) {

// Assume all configuration parameter data members have been filled in already.
void CalorimeterHitDigi::init(const dd4hep::Detector* detector, std::shared_ptr<spdlog::logger>& logger) {
m_detector = detector;
m_log = logger;

// Gaudi implements a random number generator service. It is not clear to me how this
// can work. There are multiple race conditions that occur in parallel event processing:
Expand All @@ -54,68 +48,53 @@ void CalorimeterHitDigi::AlgorithmInit(std::shared_ptr<spdlog::logger>& logger)
// now, just use default values defined in header file.

// set energy resolution numbers
m_log=logger;

if (u_eRes.empty()) {
u_eRes.resize(3);
} else if (u_eRes.size() != 3) {
m_log->error("Invalid u_eRes.size()");
throw std::runtime_error("Invalid u_eRes.size()");
if (m_cfg.eRes.empty()) {
m_cfg.eRes.resize(3);
} else if (m_cfg.eRes.size() != 3) {
m_log->error("Invalid m_cfg.eRes.size()");
throw std::runtime_error("Invalid m_cfg.eRes.size()");
}

// using juggler internal units (GeV, mm, radian, ns)
tRes = m_tRes / dd4hep::ns;
stepTDC = dd4hep::ns / m_resolutionTDC;
tRes = m_cfg.tRes / dd4hep::ns;
stepTDC = dd4hep::ns / m_cfg.resolutionTDC;

// all these are for signal sum at digitization level
merge_hits = false;
if (!u_fields.empty()) {
if (!m_cfg.fields.empty()) {
// sanity checks
if (!m_geoSvc) {
m_log->error("Unable to locate Geometry Service.");
if (!m_detector) {
m_log->error("Unable to locate geometry.");
throw std::runtime_error("Unable to locate Geometry Service.");
}
if (m_readout.empty()) {
if (m_cfg.readout.empty()) {
m_log->error("readoutClass is not provided, it is needed to know the fields in readout ids.");
throw std::runtime_error("readoutClass is not provided.");
}

// get decoders
try {
auto id_desc = m_geoSvc->detector()->readout(m_readout).idSpec();
auto id_desc = m_detector->readout(m_cfg.readout).idSpec();
id_mask = 0;
for (size_t i = 0; i < u_fields.size(); ++i) {
id_mask |= id_desc.field(u_fields[i])->mask();
for (size_t i = 0; i < m_cfg.fields.size(); ++i) {
id_mask |= id_desc.field(m_cfg.fields[i])->mask();
}
} catch (...) {
// a workaround to avoid breaking the whole analysis if a field is not in some configurations
// TODO: it should be a fatal error to not cause unexpected analysis results
m_log->warn("Failed to load ID decoder for {}, hits will not be merged.", m_readout);
// throw::runtime_error(fmt::format("Failed to load ID decoder for {}", m_readout));
m_log->warn("Failed to load ID decoder for {}, hits will not be merged.", m_cfg.readout);
// throw::runtime_error(fmt::format("Failed to load ID decoder for {}", m_cfg.readout));
return;
}
id_mask = ~id_mask;
m_log->info("ID mask in {:s}: {:#064b}", m_readout, id_mask);
m_log->info("ID mask in {:s}: {:#064b}", m_cfg.readout, id_mask);
// all checks passed
merge_hits = true;
}
}



//------------------------
// AlgorithmChangeRun
//------------------------
void CalorimeterHitDigi::AlgorithmChangeRun() {
/// This is automatically run before Process, when a new run number is seen
/// Usually we update our calibration constants by asking a JService
/// to give us the latest data for this run number
}

//------------------------
// AlgorithmProcess
//------------------------
std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::AlgorithmProcess(const edm4hep::SimCalorimeterHitCollection &simhits) {
std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::process(const edm4hep::SimCalorimeterHitCollection &simhits) {
if (merge_hits) {
return std::move(signal_sum_digi(simhits));
} else {
Expand All @@ -132,35 +111,31 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::single
const double eDep = ahit.getEnergy();

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

const double ped = m_pedMeanADC + m_normDist(generator) * m_pedSigmaADC;
const long long adc = std::llround(ped + eDep * (m_corrMeanScale + eResRel) / m_dyRangeADC * m_capADC);
const double ped = m_cfg.pedMeanADC + m_normDist(generator) * m_cfg.pedSigmaADC;
const long long adc = std::llround(ped + eDep * (m_cfg.corrMeanScale + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC);

double time = std::numeric_limits<double>::max();
for (const auto& c : ahit.getContributions()) {
if (c.getTime() <= time) {
time = c.getTime();
}
}
if (time > m_capTime) continue;
if (time > m_cfg.capTime) continue;

const long long tdc = std::llround((time + m_normDist(generator) * tRes) * stepTDC);

if (eDep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t cell ID {}", eDep, adc, time, m_capTime, tdc, ahit.getCellID());
if (eDep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {} \t cell ID {}", eDep, adc, time, m_cfg.capTime, tdc, ahit.getCellID());
rawhits->create(
ahit.getCellID(),
(adc > m_capADC ? m_capADC : adc),
(adc > m_cfg.capADC ? m_cfg.capADC : adc),
tdc
);
}
Expand All @@ -177,8 +152,8 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::signal
for (const auto &ahit : simhits) {
uint64_t hid = ahit.getCellID() & id_mask;

m_log->trace("org cell ID in {:s}: {:#064b}", m_readout, ahit.getCellID());
m_log->trace("new cell ID in {:s}: {:#064b}", m_readout, hid);
m_log->trace("org cell ID in {:s}: {:#064b}", m_cfg.readout, ahit.getCellID());
m_log->trace("new cell ID in {:s}: {:#064b}", m_cfg.readout, hid);

merge_map[hid].push_back(ix);

Expand All @@ -202,7 +177,7 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::signal
timeC = c.getTime();
}
}
if (timeC > m_capTime) continue;
if (timeC > m_cfg.capTime) continue;
edep += hit.getEnergy();
m_log->trace("adding {} \t total: {}", hit.getEnergy(), edep);

Expand All @@ -221,29 +196,25 @@ std::unique_ptr<edm4hep::RawCalorimeterHitCollection> CalorimeterHitDigi::signal
}
}

// double eResRel = 0.;
// safety check
const double eResRel = (edep > m_threshold)
? m_normDist(generator) * u_eRes[0] / std::sqrt(edep) +
m_normDist(generator) * u_eRes[1] +
m_normDist(generator) * u_eRes[2] / edep
const double eResRel = (edep > m_cfg.threshold)
? m_normDist(generator) * m_cfg.eRes[0] / std::sqrt(edep) +
m_normDist(generator) * m_cfg.eRes[1] +
m_normDist(generator) * m_cfg.eRes[2] / edep
: 0;
// if (edep > 1e-6) {
// eResRel = m_normDist(generator) * u_eRes[0] / std::sqrt(edep) +
// m_normDist(generator) * u_eRes[1] +
// m_normDist(generator) * u_eRes[2] / edep;
// }
double ped = m_pedMeanADC + m_normDist(generator) * m_pedSigmaADC;
unsigned long long adc = std::llround(ped + edep * (m_corrMeanScale + eResRel) / m_dyRangeADC * m_capADC);
double ped = m_cfg.pedMeanADC + m_normDist(generator) * m_cfg.pedSigmaADC;
unsigned long long adc = std::llround(ped + edep * (m_cfg.corrMeanScale + eResRel) / m_cfg.dyRangeADC * m_cfg.capADC);
unsigned long long tdc = std::llround((time + m_normDist(generator) * tRes) * stepTDC);

if (edep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {}", edep, adc, time, m_capTime, tdc);
if (edep> 1.e-3) m_log->trace("E sim {} \t adc: {} \t time: {}\t maxtime: {} \t tdc: {}", edep, adc, time, m_cfg.capTime, tdc);
rawhits->create(
mid,
(adc > m_capADC ? m_capADC : adc),
(adc > m_cfg.capADC ? m_cfg.capADC : adc),
tdc
);
}

return std::move(rawhits);
}

} // namespace eicrecon
66 changes: 18 additions & 48 deletions src/algorithms/calorimetry/CalorimeterHitDigi.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,72 +16,42 @@
#include <memory>
#include <random>

#include "services/geometry/dd4hep/JDD4hep_service.h"
#include <DD4hep/Detector.h>

#include <edm4hep/SimCalorimeterHitCollection.h>
#include <edm4hep/RawCalorimeterHitCollection.h>
#include <spdlog/spdlog.h>

class CalorimeterHitDigi {
#include "algorithms/interfaces/WithPodConfig.h"
#include "CalorimeterHitDigiConfig.h"

// Insert any member variables here
namespace eicrecon {

public:
CalorimeterHitDigi() = default;
void AlgorithmInit(std::shared_ptr<spdlog::logger>& logger);
void AlgorithmChangeRun() ;
std::unique_ptr<edm4hep::RawCalorimeterHitCollection> AlgorithmProcess(const edm4hep::SimCalorimeterHitCollection &simhits) ;
class CalorimeterHitDigi : public WithPodConfig<CalorimeterHitDigiConfig> {

//-------- Configuration Parameters ------------
//instantiate new spdlog logger
std::shared_ptr<spdlog::logger> m_log;
public:
void init(const dd4hep::Detector* detector, std::shared_ptr<spdlog::logger>& logger);
std::unique_ptr<edm4hep::RawCalorimeterHitCollection> process(const edm4hep::SimCalorimeterHitCollection &simhits) ;

// Name of input data type (collection)
std::string m_input_tag;

// additional smearing resolutions
std::vector<double> u_eRes;
double m_tRes;

// single hit energy deposition threshold
double m_threshold=1.0*dd4hep::keV; // {this, "threshold", 1. * keV};

// digitization settings
unsigned int m_capADC;
double m_capTime = 1000.;
double m_dyRangeADC;
unsigned int m_pedMeanADC;
double m_pedSigmaADC;
double m_resolutionTDC;
double m_corrMeanScale;

// signal sums
std::vector<std::string> u_fields;
std::string m_geoSvcName;
std::string m_readout;

// This may be used to declare the data members as JANA configuration parameters.
// This should compile OK even without JANA so long as you don't try using it.
// To use it, do something like the following:
//
// mycalohitdigi->SetJANAConfigParameters( japp, "BEMC");
//
// The above will register config. parameters like: "BEMC:tag".
// The configuration parameter members of this class should be set to their
// defaults *before* calling this.
//-----------------------------------------------
private:

// unitless counterparts of inputs
double dyRangeADC{0}, stepTDC{0}, tRes{0};
// variables for merging at digitization step
bool merge_hits = false;
std::shared_ptr<JDD4hep_service> m_geoSvc;

uint64_t id_mask{0};

private:
private:
const dd4hep::Detector* m_detector;
std::shared_ptr<spdlog::logger> m_log;

std::default_random_engine generator; // TODO: need something more appropriate here
std::normal_distribution<double> m_normDist; // defaults to mean=0, sigma=1

std::unique_ptr<edm4hep::RawCalorimeterHitCollection> single_hits_digi(const edm4hep::SimCalorimeterHitCollection &simhits);
std::unique_ptr<edm4hep::RawCalorimeterHitCollection> signal_sum_digi(const edm4hep::SimCalorimeterHitCollection &simhits);
};

};

} // namespace eicrecon
34 changes: 34 additions & 0 deletions src/algorithms/calorimetry/CalorimeterHitDigiConfig.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,34 @@
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2023 Wouter Deconinck

#pragma once

#include <string>
#include <vector>

namespace eicrecon {

struct CalorimeterHitDigiConfig {

std::vector<double> eRes;
double tRes;

// single hit energy deposition threshold
double threshold{1.0*dd4hep::keV};

// digitization settings
unsigned int capADC{1};
double capTime{1000}; // dynamic range in ns
double dyRangeADC{1};
unsigned int pedMeanADC{0};
double pedSigmaADC{0};
double resolutionTDC{1};
double corrMeanScale{1};

// signal sums
std::string readout{""};
std::vector<std::string> fields{};

};

} // eicrecon
Loading

0 comments on commit 20c9d80

Please sign in to comment.