Skip to content

Commit

Permalink
feat: use digitization algorithms from EICrecon
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdconinc
wdconinc committed Jul 27, 2024
1 parent dd4ec49 commit 3f59618
Show file tree
Hide file tree
Showing 3 changed files with 10 additions and 317 deletions.
2 changes: 2 additions & 0 deletions JugDigi/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,12 @@ gaudi_add_module(JugDigiPlugins
JugBase
JugAlgo
EICrecon::algorithms_calorimetry_library
EICrecon::algorithms_digi_library
ROOT::Core ROOT::RIO ROOT::Tree
EDM4HEP::edm4hep
EDM4EIC::edm4eic
k4FWCore::k4FWCore
GENCONF_OPTIONS --load-library=JugAlgoPreload.so
)

target_include_directories(JugDigiPlugins PUBLIC
Expand Down
221 changes: 4 additions & 217 deletions JugDigi/src/components/PhotoMultiplierDigi.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,221 +1,8 @@
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2022 Chao Peng
// Copyright (C) 2024 Wouter Deconinck

/* General PhotoMultiplier Digitization
*
* Apply the given quantum efficiency for photon detection
* Converts the number of detected photons to signal amplitude
*
* Author: Chao Peng (ANL)
* Date: 10/02/2020
*/

#include <iterator>
#include <algorithm>
#include <unordered_map>
#include <cmath>

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

#include <k4FWCore/DataHandle.h>

// Event Model related classes
#include "edm4eic/RawTrackerHitCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"


using namespace Gaudi::Units;

namespace Jug::Digi {

/** PhotoMultiplierDigi.
*
* \ingroup digi
*/
class PhotoMultiplierDigi : public GaudiAlgorithm
{
public:
DataHandle<edm4hep::SimTrackerHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<edm4eic::RawTrackerHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
Gaudi::Property<std::vector<std::pair<double, double>>>
u_quantumEfficiency{this, "quantumEfficiency", {{2.6*eV, 0.3}, {7.0*eV, 0.3}}};
Gaudi::Property<double> m_hitTimeWindow{this, "hitTimeWindow", 20.0*ns};
Gaudi::Property<double> m_timeStep{this, "timeStep", 0.0625*ns};
Gaudi::Property<double> m_speMean{this, "speMean", 80.0};
Gaudi::Property<double> m_speError{this, "speError", 16.0};
Gaudi::Property<double> m_pedMean{this, "pedMean", 200.0};
Gaudi::Property<double> m_pedError{this, "pedError", 3.0};
Rndm::Numbers m_rngUni, m_rngNorm;

// constructor
PhotoMultiplierDigi(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;
}

auto randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
auto sc1 = m_rngUni.initialize(randSvc, Rndm::Flat(0., 1.));
auto sc2 = m_rngNorm.initialize(randSvc, Rndm::Gauss(0., 1.));
if (!sc1.isSuccess() || !sc2.isSuccess()) {
error() << "Cannot initialize random generator!" << endmsg;
return StatusCode::FAILURE;
}

qe_init();

return StatusCode::SUCCESS;
}

StatusCode execute() override
{
// input collection
const auto &sim = *m_inputHitCollection.get();
// Create output collections
auto &raw = *m_outputHitCollection.createAndPut();

struct HitData { int npe; double signal; double time; };
std::unordered_map<decltype(edm4eic::RawTrackerHitData::cellID), std::vector<HitData>> hit_groups;
// collect the photon hit in the same cell
// calculate signal
for(const auto& ahit : sim) {
// quantum efficiency
if (!qe_pass(ahit.getEDep(), m_rngUni())) {
continue;
}
// cell id, time, signal amplitude
uint64_t id = ahit.getCellID();
double time = ahit.getMCParticle().getTime();
double amp = m_speMean + m_rngNorm()*m_speError;

// group hits
auto it = hit_groups.find(id);
if (it != hit_groups.end()) {
size_t i = 0;
for (auto git = it->second.begin(); git != it->second.end(); ++git, ++i) {
if (std::abs(time - git->time) <= (m_hitTimeWindow/ns)) {
git->npe += 1;
git->signal += amp;
break;
}
}
// no hits group found
if (i >= it->second.size()) {
it->second.emplace_back(HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time});
}
} else {
hit_groups[id] = {HitData{1, amp + m_pedMean + m_pedError*m_rngNorm(), time}};
}
}

// build hit
for (auto &it : hit_groups) {
for (auto &data : it.second) {
raw.create(
it.first,
static_cast<decltype(edm4eic::RawTrackerHitData::charge)>(data.signal),
static_cast<decltype(edm4eic::RawTrackerHitData::timeStamp)>(data.time/(m_timeStep/ns))
);
}
}

return StatusCode::SUCCESS;
}

private:
void qe_init()
{
auto &qeff = u_quantumEfficiency.value();

// sort quantum efficiency data first
std::sort(qeff.begin(), qeff.end(),
[] (const std::pair<double, double> &v1, const std::pair<double, double> &v2) {
return v1.first < v2.first;
});

// sanity checks
if (qeff.empty()) {
qeff = {{2.6*eV, 0.3}, {7.0*eV, 0.3}};
warning() << "Invalid quantum efficiency data provided, using default values: " << qeff << endmsg;
}
if (qeff.front().first > 3.0*eV) {
warning() << "Quantum efficiency data start from " << qeff.front().first/eV
<< " eV, maybe you are using wrong units?" << endmsg;
}
if (qeff.back().first < 6.0*eV) {
warning() << "Quantum efficiency data end at " << qeff.back().first/eV
<< " eV, maybe you are using wrong units?" << endmsg;
}
}

// helper function for linear interpolation
// Comp return is defined as: equal, 0; greater, > 0; less, < 0
template<class RndmIter, typename T, class Compare>
RndmIter interval_search(RndmIter beg, RndmIter end, const T &val, Compare comp) const
{
// special cases
auto dist = std::distance(beg, end);
if ((dist < 2) || (comp(*beg, val) > 0) || (comp(*std::prev(end), val) < 0)) {
return end;
}
auto mid = std::next(beg, dist / 2);

while (mid != end) {
if (comp(*mid, val) == 0) {
return mid;
} else if (comp(*mid, val) > 0) {
end = mid;
} else {
beg = std::next(mid);
}
mid = std::next(beg, std::distance(beg, end)/2);
}

if (mid == end || comp(*mid, val) > 0) {
return std::prev(mid);
}
return mid;
}

bool qe_pass(double ev, double rand) const
{
const auto &qeff = u_quantumEfficiency.value();
auto it = interval_search(qeff.begin(), qeff.end(), ev,
[] (const std::pair<double, double> &vals, double val) {
return vals.first - val;
});

if (it == qeff.end()) {
// info() << ev/eV << " eV is out of QE data range, assuming 0% efficiency" << endmsg;
return false;
}

double prob = it->second;
auto itn = std::next(it);
if (itn != qeff.end() && (itn->first - it->first != 0)) {
prob = (it->second*(itn->first - ev) + itn->second*(ev - it->first)) / (itn->first - it->first);
}

// info() << ev/eV << " eV, QE: " << prob*100. << "%" << endmsg;
return rand <= prob;
}
};
#include <JugAlgo/Algorithm.h>
#include <EICrecon/algorithms/digi/PhotoMultiplierHitDigi.h>

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

} // namespace Jug::Digi
JUGALGO_DEFINE_ALGORITHM(PhotoMultiplierDigi, eicrecon::PhotoMultiplierHitDigi, Jug::Digi)
104 changes: 4 additions & 100 deletions JugDigi/src/components/SiliconTrackerDigi.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,104 +1,8 @@
// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Sylvester Joosten
// Copyright (C) 2024 Wouter Deconinck

#include <algorithm>
#include <cmath>
#include <JugAlgo/Algorithm.h>
#include <EICrecon/algorithms/digi/SiliconTrackerDigi.h>

#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 <k4FWCore/DataHandle.h>

// Event Model related classes
// edm4hep's tracker hit is the input collectiopn
#include "edm4hep/MCParticle.h"
#include "edm4hep/SimTrackerHitCollection.h"
// edm4eic's RawTrackerHit is the output
#include "edm4eic/RawTrackerHitCollection.h"

namespace Jug::Digi {

/** Silicon detector digitization.
*
* \ingroup digi
*/
class SiliconTrackerDigi : public GaudiAlgorithm {
private:
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10}; // todo : add units
Gaudi::Property<double> m_threshold{this, "threshold", 0. * Gaudi::Units::keV};
Rndm::Numbers m_gaussDist;
DataHandle<edm4hep::SimTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<edm4eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};

public:
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
SiliconTrackerDigi(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;
}
IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(0.0, m_timeResolution.value()));
if (!sc.isSuccess()) {
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override {
// input collection
const auto* const simhits = m_inputHitCollection.get();
// Create output collections
auto* rawhits = m_outputHitCollection.createAndPut();
// edm4eic::RawTrackerHitCollection* rawHitCollection = new edm4eic::RawTrackerHitCollection();
std::map<long long, int> cell_hit_map;
for (const auto& ahit : *simhits) {
if (msgLevel(MSG::DEBUG)) {
debug() << "--------------------" << ahit.getCellID() << endmsg;
debug() << "Hit in cellID = " << ahit.getCellID() << endmsg;
debug() << " position = (" << ahit.getPosition().x << "," << ahit.getPosition().y << ","
<< ahit.getPosition().z << ")" << endmsg;
debug() << " xy_radius = " << std::hypot(ahit.getPosition().x, ahit.getPosition().y) << endmsg;
debug() << " momentum = (" << ahit.getMomentum().x << "," << ahit.getMomentum().y << ","
<< ahit.getMomentum().z << ")" << endmsg;
}
if (ahit.getEDep() * Gaudi::Units::keV < m_threshold) {
if (msgLevel(MSG::DEBUG)) {
debug() << " edep = " << ahit.getEDep() << " (below threshold of " << m_threshold / Gaudi::Units::keV
<< " keV)" << endmsg;
}
continue;
} else {
if (msgLevel(MSG::DEBUG)) {
debug() << " edep = " << ahit.getEDep() << endmsg;
}
}
if (cell_hit_map.count(ahit.getCellID()) == 0) {
cell_hit_map[ahit.getCellID()] = rawhits->size();
rawhits->create(
ahit.getCellID(),
ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3, // ns->fs
std::llround(ahit.getEDep() * 1e6)
);
} else {
auto hit = (*rawhits)[cell_hit_map[ahit.getCellID()]];
hit.setTimeStamp(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
auto ch = hit.getCharge();
hit.setCharge(ch + std::llround(ahit.getEDep() * 1e6));
}
}
return StatusCode::SUCCESS;
}
};
// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
DECLARE_COMPONENT(SiliconTrackerDigi)

} // namespace Jug::Digi
JUGALGO_DEFINE_ALGORITHM(SiliconTrackerDigi, eicrecon::SiliconTrackerDigi, Jug::Digi)

0 comments on commit 3f59618

Please sign in to comment.