Skip to content

GitLab

Commit 2e4368e2 authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

Support pure edm4hep output from ddsim

parent 25012406
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
View file @ 2e4368e2
......@@ -68,13 +68,7 @@ juggler:local:
image: eicweb.phy.anl.gov:4567/containers/eic_container/jug_xl:nightly
stage: build
script:
## first install NPDET and EICD to ensure the latest version, then build juggler
- |
git clone https://eicweb.phy.anl.gov/eic/npdet.git /tmp/npdet
cd /tmp/npdet && git checkout $JUGGLER_NPDET_VERSION && cd -
cmake -B /tmp/build -S /tmp/npdet -DCMAKE_INSTALL_PREFIX=/usr/local
cmake --build /tmp/build -j40 -- install
rm -rf /tmp/build /tmp/npdet
## first install EICD to ensure the latest version, then build juggler
- |
git clone https://eicweb.phy.anl.gov/eic/eicd.git /tmp/eicd
cd /tmp/eicd && git checkout $JUGGLER_EICD_VERSION && cd -
......
CMakeLists.txt
View file @ 2e4368e2
......@@ -23,7 +23,7 @@ if(ENABLE_CLANG_TIDY)
endif()
find_package(EICD REQUIRED)
find_package(NPDet REQUIRED)
find_package(EDM4HEP REQUIRED)
find_package(podio 0.11...0.14 REQUIRED)
add_definitions("-Dpodio_VERSION_MAJOR=${podio_VERSION_MAJOR}")
......
JugBase/CMakeLists.txt
View file @ 2e4368e2
......@@ -42,7 +42,7 @@ gaudi_add_module(JugBasePlugins
Gaudi::GaudiKernel Gaudi::GaudiAlgLib
ROOT::Core ROOT::RIO ROOT::Tree
JugBase
NPDet::DD4podIO
EDM4HEP::edm4hep
DD4hep::DDRec
ActsCore ActsPluginDD4hep ActsPluginJson
EICD::eicd
......
JugBase/src/components/InputCopier.cpp
View file @ 2e4368e2
......@@ -10,10 +10,9 @@
#include "JugBase/DataHandle.h"
// Event Model related classes
//#include "GaudiExamples/MyTrack.h"
#include "dd4pod/Geant4ParticleCollection.h"
#include "dd4pod/CalorimeterHitCollection.h"
#include "dd4pod/TrackerHitCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"
namespace Jug {
......@@ -29,7 +28,7 @@ namespace Jug {
public:
InputCopier(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputCollection", m_inputHitCollection, "mcparticles");
declareProperty("inputCollection", m_inputHitCollection, "MCParticles");
declareProperty("outputCollection", m_outputHitCollection, "genparticles");
}
StatusCode initialize() override
......@@ -55,58 +54,19 @@ namespace Jug {
return StatusCode::SUCCESS;
}
DataHandle<T_IN> m_inputHitCollection{"mcparticles", Gaudi::DataHandle::Reader, this};
DataHandle<T_IN> m_inputHitCollection{"MCParticles", Gaudi::DataHandle::Reader, this};
DataHandle<T_OUT> m_outputHitCollection{"genparticles", Gaudi::DataHandle::Writer, this};
};
using CalorimeterColCopier = InputCopier<dd4pod::CalorimeterHitCollection, dd4pod::CalorimeterHitCollection>;
using CalorimeterColCopier = InputCopier<edm4hep::SimCalorimeterHitCollection, edm4hep::SimCalorimeterHitCollection>;
DECLARE_COMPONENT(CalorimeterColCopier)
using TrackerColCopier = InputCopier<dd4pod::TrackerHitCollection, dd4pod::TrackerHitCollection>;
using TrackerColCopier = InputCopier<edm4hep::SimTrackerHitCollection, edm4hep::SimTrackerHitCollection>;
DECLARE_COMPONENT(TrackerColCopier)
using MCCopier = InputCopier<dd4pod::Geant4ParticleCollection, dd4pod::Geant4ParticleCollection>;
using MCCopier = InputCopier<edm4hep::MCParticleCollection, edm4hep::MCParticleCollection>;
DECLARE_COMPONENT(MCCopier)
//class MCCopier : public GaudiAlgorithm {
//public:
// // ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
// MCCopier(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
// {
// declareProperty("inputCollection", m_inputHitCollection, "mcparticles");
// declareProperty("outputCollection", m_outputHitCollection, "genparticles");
// }
// StatusCode initialize() override
// {
// if (GaudiAlgorithm::initialize().isFailure())
// return StatusCode::FAILURE;
// // f_counter = m_starting_value.value();
// return StatusCode::SUCCESS;
// }
// StatusCode execute() override
// {
// // input collection
// const dd4pod::Geant4ParticleCollection* simhits = m_inputHitCollection.get();
// // Create output collections
// auto out_parts = m_outputHitCollection.createAndPut();
// // std::copy(std::begin(*simhits),std::end(*simhits),std::begin(*out_parts));
// for (const auto& ahit : *simhits) {
// // std::cout << ahit << "\n";
// // eic::RawCalorimeterHit rawhit((long long)ahit.cellID(),
// // (long long)ahit.energyDeposit() * 100, 0);
// // dd4pod::Geant4Particle bpart = ahit;// dd4pod::Geant4ParticleCollection():
// out_parts->push_back(ahit.clone());
// }
// return StatusCode::SUCCESS;
// }
// DataHandle<dd4pod::Geant4ParticleCollection> m_inputHitCollection{"mcparticles", Gaudi::DataHandle::Reader, this};
// DataHandle<dd4pod::Geant4ParticleCollection> m_outputHitCollection{"genparticles", Gaudi::DataHandle::Writer,
// this};
//};
//DECLARE_COMPONENT(MCCopier)
} // namespace Examples
} // namespace Gaudi
JugBase/src/components/MC2DummyParticle.cpp
View file @ 2e4368e2
......@@ -10,26 +10,26 @@
#include "JugBase/UniqueID.h"
// Event Model related classes
#include "dd4pod/Geant4ParticleCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
namespace Jug::Base {
class MC2DummyParticle : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
public:
DataHandle<dd4pod::Geant4ParticleCollection> m_inputParticles{"mcparticles", Gaudi::DataHandle::Reader, this};
DataHandle<edm4hep::MCParticleCollection> m_inputParticles{"MCParticles", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ReconstructedParticleCollection> m_outputParticles{"DummyReconstructedParticles",
Gaudi::DataHandle::Writer, this};
Rndm::Numbers m_gaussDist;
Gaudi::Property<double> m_smearing{this, "smearing", 0.01 /* 1 percent*/};
const int32_t kMonteCarloSource{uniqueID<int32_t>("mcparticles")};
const int32_t kMonteCarloSource{uniqueID<int32_t>("MCParticles")};
MC2DummyParticle(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
, AlgorithmIDMixin(name, info())
{
declareProperty("inputCollection", m_inputParticles, "mcparticles");
declareProperty("inputCollection", m_inputParticles, "MCParticles");
declareProperty("outputCollection", m_outputParticles, "DummyReconstructedParticles");
}
StatusCode initialize() override
......@@ -47,14 +47,14 @@ namespace Jug::Base {
StatusCode execute() override
{
// input collection
const dd4pod::Geant4ParticleCollection* parts = m_inputParticles.get();
const edm4hep::MCParticleCollection* parts = m_inputParticles.get();
// output collection
auto& out_parts = *(m_outputParticles.createAndPut());
int ID = 0;
for (const auto& p : *parts) {
if (p.genStatus() > 1) {
if (p.getGeneratorStatus() > 1) {
if (msgLevel(MSG::DEBUG)) {
debug() << "ignoring particle with genStatus = " << p.genStatus() << endmsg;
debug() << "ignoring particle with generatorStatus = " << p.getGeneratorStatus() << endmsg;
}
continue;
}
......@@ -62,33 +62,34 @@ namespace Jug::Base {
// for now just use total momentum smearing as this is the largest effect,
// ideally we should also smear the angles but this should be good enough
// for now.
const double pgen = p.ps().mag();
const float momentum = pgen * m_gaussDist();
const float energy = p.energy();
const float px = p.ps().x * momentum / pgen;
const float py = p.ps().y * momentum / pgen;
const float pz = p.ps().z * momentum / pgen;
const auto pvec = p.getMomentum();
const auto pgen = std::hypot(pvec.x, pvec.y, pvec.z);
const auto momentum = pgen * m_gaussDist();
const auto energy = p.getEnergy();
const auto px = p.getMomentum().x * momentum / pgen;
const auto py = p.getMomentum().y * momentum / pgen;
const auto pz = p.getMomentum().z * momentum / pgen;
// @TODO: vertex smearing
const float vx = p.vs().x;
const float vy = p.vs().y;
const float vz = p.vs().z;
const auto vx = p.getVertex().x;
const auto vy = p.getVertex().y;
const auto vz = p.getVertex().z;
eic::VectorXYZ psmear{px, py, pz};
const auto p_phi = std::atan2(py, px);
const auto p_theta = std::atan2(std::hypot(px, py), pz);
auto rec_part = out_parts.create();
rec_part.ID({ID++, algorithmID()});
rec_part.p(psmear);
rec_part.p({px, py, pz});
rec_part.v({vx, vy, vz});
rec_part.time(static_cast<float>(p.time()));
rec_part.pid(p.pdgID());
rec_part.status(p.genStatus());
rec_part.charge(p.charge());
rec_part.time(p.getTime());
rec_part.pid(p.getPDG());
rec_part.status(p.getGeneratorStatus());
rec_part.charge(p.getCharge());
rec_part.weight(1.);
rec_part.direction({psmear.theta(), psmear.phi()});
rec_part.direction({p_theta, p_phi});
rec_part.momentum(momentum);
rec_part.energy(energy);
rec_part.mass(p.mass());
rec_part.mcID({p.ID(), kMonteCarloSource});
rec_part.mass(p.getMass());
}
return StatusCode::SUCCESS;
}
......
JugBase/src/components/ReadTestConsumer.cxx
View file @ 2e4368e2
......@@ -2,13 +2,13 @@
#include "JugBase/DataHandle.h"
#include "dd4pod/Geant4ParticleCollection.h"
#include "edm4hep/MCParticleCollection.h"
class ReadTestConsumer : public GaudiAlgorithm {
public:
ReadTestConsumer(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc), m_genParticles("mcparticles", Gaudi::DataHandle::Reader, this) {
: GaudiAlgorithm(name, svcLoc), m_genParticles("MCParticles", Gaudi::DataHandle::Reader, this) {
declareProperty("genParticles", m_genParticles, "mc particles to read");
}
......@@ -20,7 +20,7 @@ public:
StatusCode execute() {
// Read the input
const dd4pod::Geant4ParticleCollection* mcparticles = m_genParticles.get();
const edm4hep::MCParticleCollection* mcparticles = m_genParticles.get();
// Does the reading work?
debug() << mcparticles << endmsg;
......@@ -40,6 +40,6 @@ public:
private:
/// Particles to read
DataHandle<dd4pod::Geant4ParticleCollection> m_genParticles;
DataHandle<edm4hep::MCParticleCollection> m_genParticles;
};
DECLARE_COMPONENT(ReadTestConsumer)
JugBase/tests/options/col_copier.py
View file @ 2e4368e2
......@@ -3,7 +3,7 @@ from Configurables import ApplicationMgr, EICDataSvc, PodioOutput
podioevent = EICDataSvc("EventDataSvc", inputs=["derp.root"], OutputLevel=DEBUG)
from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
from Configurables import Jug__Base__InputCopier_edm4hep__MCParticleCollection_edm4hep__MCParticleCollection_ as MCCopier
from Configurables import PodioInput, ReadTestConsumer
podioinput = PodioInput("PodioReader", collections=["mcparticles"], OutputLevel=DEBUG)
checker = ReadTestConsumer()
......
JugDigi/CMakeLists.txt
View file @ 2e4368e2
......@@ -14,7 +14,7 @@ gaudi_add_module(JugDigiPlugins
Gaudi::GaudiAlgLib Gaudi::GaudiKernel
JugBase
ROOT::Core ROOT::RIO ROOT::Tree
NPDet::DD4podIO
EDM4HEP::edm4hep
EICD::eicd
)
......
JugDigi/src/components/CalorimeterBirksCorr.cpp
View file @ 2e4368e2
// Apply Birks Law to correct the energy deposit
// It uses the contributions member in dd4pod::CalorimeterHit, so simulation must enable storeCalorimeterContributions
// It uses the contributions member in edm4hep::CalorimeterHit, so simulation must enable storeCalorimeterContributions
//
// Author: Chao Peng
// Date: 09/29/2021
......@@ -18,7 +18,7 @@
#include "JugBase/IParticleSvc.h"
// Event Model related classes
#include "dd4pod/CalorimeterHitCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitData.h"
......@@ -38,9 +38,9 @@ namespace Jug::Digi {
// digitization settings
Gaudi::Property<double> m_birksConstant{this, "birksConstant", 0.126*mm/MeV};
DataHandle<dd4pod::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
DataHandle<edm4hep::SimCalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<dd4pod::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
DataHandle<edm4hep::SimCalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};
SmartIF<IParticleSvc> m_pidSvc;
......@@ -80,22 +80,21 @@ namespace Jug::Digi {
{
auto& ohits = *m_outputHitCollection.createAndPut();
for (const auto& hit : *m_inputHitCollection.get()) {
double lightyield = 0.;
for (auto &truth : hit.contributions()) {
const double charge = m_pidSvc->particle(truth.pdgID).charge;
auto ohit = ohits->create(hit.getCellID(), hit.getEnergy(), hit.getPosition());
double energy = 0.;
for (const auto &c: hit.getContributions()) {
ohit.addToContributions(c);
const double charge = m_pidSvc->particle(c.getPDG()).charge;
// some tolerance for precision
if (std::abs(charge) > 1e-5) {
lightyield += truth.deposit / (1. + truth.deposit / truth.length * birksConstant);
// FIXME
//energy += c.getEnergy() / (1. + c.getEnergy() / c.length * birksConstant);
error() << "edm4hep::CaloHitContribution has no length field for Birks correction." << endmsg;
return StatusCode::FAILURE;
}
}
auto ohit = ohits->create();
ohit.cellID(hit.cellID());
ohit.flag(hit.flag());
ohit.g4ID(hit.g4ID());
ohit.position(hit.position());
ohit.truth(hit.truth());
// replace energy deposit with Birks Law corrected value
ohit.energyDeposit(lightyield);
ohit.setEnergy(energy);
}
return StatusCode::SUCCESS;
}
......
JugDigi/src/components/CalorimeterHitDigi.cpp
View file @ 2e4368e2
......@@ -28,7 +28,7 @@
#include "fmt/ranges.h"
// Event Model related classes
#include "dd4pod/CalorimeterHitCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitData.h"
......@@ -76,7 +76,7 @@ namespace Jug::Digi {
SmartIF<IGeoSvc> m_geoSvc;
uint64_t id_mask, ref_mask;
DataHandle<dd4pod::CalorimeterHitCollection> m_inputHitCollection{
DataHandle<edm4hep::SimCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawCalorimeterHitCollection> m_outputHitCollection{
"outputHitCollection", Gaudi::DataHandle::Writer, this};
......@@ -171,7 +171,7 @@ namespace Jug::Digi {
int nhits = 0;
for (const auto& ahit : *simhits) {
// Note: juggler internal unit of energy is GeV
const double eDep = ahit.energyDeposit();
const double eDep = ahit.getEnergy();
// apply additional calorimeter noise to corrected energy deposit
const double eResRel = (eDep > 1e-6)
......@@ -181,11 +181,16 @@ namespace Jug::Digi {
const double ped = m_pedMeanADC + m_normDist() * m_pedSigmaADC;
const long long adc = std::llround(ped + m_corrMeanScale * eDep * (1. + eResRel) / dyRangeADC * m_capADC);
const long long tdc = std::llround((ahit.truth().time + m_normDist() * tRes) * stepTDC);
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);
eic::RawCalorimeterHit rawhit(
{nhits++, algorithmID()},
(long long)ahit.cellID(),
(long long)ahit.getCellID(),
(adc > m_capADC.value() ? m_capADC.value() : adc),
tdc
);
......@@ -198,9 +203,9 @@ namespace Jug::Digi {
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;
std::unordered_map<long long, std::vector<edm4hep::ConstSimCalorimeterHit>> merge_map;
for (const auto &ahit : *simhits) {
int64_t hid = (ahit.cellID() & id_mask) | ref_mask;
int64_t hid = (ahit.getCellID() & id_mask) | ref_mask;
auto it = merge_map.find(hid);
if (it == merge_map.end()) {
......@@ -213,15 +218,17 @@ namespace Jug::Digi {
// 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();
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].energyDeposit();
if (hits[i].energyDeposit() > max_edep) {
max_edep = hits[i].energyDeposit();
time = hits[i].truth().time;
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();
}
}
......
JugDigi/src/components/PhotoMultiplierDigi.cpp
View file @ 2e4368e2
......@@ -22,7 +22,8 @@
// Event Model related classes
#include "eicd/RawPMTHitCollection.h"
#include "dd4pod/TrackerHitCollection.h"
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimTrackerHitCollection.h"
using namespace Gaudi::Units;
......@@ -36,7 +37,7 @@ namespace Jug::Digi {
class PhotoMultiplierDigi : public GaudiAlgorithm, AlgorithmIDMixin<>
{
public:
DataHandle<dd4pod::TrackerHitCollection>
DataHandle<edm4hep::SimTrackerHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawPMTHitCollection>
m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
......@@ -91,12 +92,12 @@ public:
// calculate signal
for(const auto& ahit : sim) {
// quantum efficiency
if (!qe_pass(ahit.energyDeposit(), m_rngUni())) {
if (!qe_pass(ahit.getEDep(), m_rngUni())) {
continue;
}
// cell id, time, signal amplitude
long long id = ahit.cellID();
double time = ahit.truth().time;
long long id = ahit.getCellID();
double time = ahit.getMCParticle().getTime();
double amp = m_speMean + m_rngNorm()*m_speError;
// group hits
......
JugDigi/src/components/SiliconTrackerDigi.cpp
View file @ 2e4368e2
......@@ -11,8 +11,9 @@
#include "JugBase/UniqueID.h"
// Event Model related classes
// dd4pod's tracker hit is the input collectiopn
#include "dd4pod/TrackerHitCollection.h"
// edm4hep's tracker hit is the input collectiopn
#include "edm4hep/MCParticle.h"
#include "edm4hep/SimTrackerHitCollection.h"
// eicd's RawTrackerHit is the output
#include "eicd/RawTrackerHitCollection.h"
......@@ -27,7 +28,7 @@ namespace Jug::Digi {
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<dd4pod::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
DataHandle<edm4hep::SimTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};
......@@ -55,7 +56,7 @@ namespace Jug::Digi {
StatusCode execute() override
{
// input collection
const dd4pod::TrackerHitCollection* simhits = m_inputHitCollection.get();
auto simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
// eic::RawTrackerHitCollection* rawHitCollection = new eic::RawTrackerHitCollection();
......@@ -63,35 +64,34 @@ namespace Jug::Digi {
int ID = 0;
for (const auto& ahit : *simhits) {
if (msgLevel(MSG::DEBUG)) {
debug() << "--------------------" << ahit.cellID() << endmsg;
debug() << "Hit in cellID = " << ahit.cellID() << endmsg;
debug() << " position = (" << ahit.position().x << "," << ahit.position().y <<","<< ahit.position().z << ")" << endmsg;
debug() << " xy_radius = " << std::hypot(ahit.position().x , ahit.position().y ) << endmsg;
debug() << " momentum = (" << ahit.momentum().x << "," << ahit.momentum().y <<","<< ahit.momentum().z << ")" << endmsg;
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.energyDeposit() * Gaudi::Units::keV < m_threshold) {
if (ahit.getEDep() * Gaudi::Units::keV < m_threshold) {
if (msgLevel(MSG::DEBUG)) {
debug() << " edep = " << ahit.energyDeposit() << " (below threshold of " << m_threshold / Gaudi::Units::keV << " keV)" << endmsg;
debug() << " edep = " << ahit.getEDep() << " (below threshold of " << m_threshold / Gaudi::Units::keV << " keV)" << endmsg;
}
continue;
} else {
if (msgLevel(MSG::DEBUG)) {
debug() << " edep = " << ahit.energyDeposit() << endmsg;
debug() << " edep = " << ahit.getEDep() << endmsg;
}
}
// std::cout << ahit << "\n";
if (cell_hit_map.count(ahit.cellID()) == 0) {
cell_hit_map[ahit.cellID()] = rawhits->size();
if (cell_hit_map.count(ahit.getCellID()) == 0) {
cell_hit_map[ahit.getCellID()] = rawhits->size();
eic::RawTrackerHit rawhit({ID++, algorithmID()},
(int64_t)ahit.cellID(),
ahit.truth().time * 1e6 + m_gaussDist() * 1e3, // ns->fs
std::llround(ahit.energyDeposit() * 1e6));
(int64_t)ahit.getCellID(),
ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3, // ns->fs
std::llround(ahit.getEDep() * 1e6));
rawhits->push_back(rawhit);
} else {
auto hit = (*rawhits)[cell_hit_map[ahit.cellID()]];
hit.time(ahit.truth().time * 1e6 + m_gaussDist() * 1e3);
auto hit = (*rawhits)[cell_hit_map[ahit.getCellID()]];
hit.time(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
auto ch = hit.charge();
hit.charge(ch + std::llround(ahit.energyDeposit() * 1e6));
hit.charge(ch + std::llround(ahit.getEDep() * 1e6));
}
}
return StatusCode::SUCCESS;
......
JugDigi/src/components/UFSDTrackerDigi.cpp
View file @ 2e4368e2
......@@ -11,10 +11,9 @@
#include "JugBase/UniqueID.h"
// Event Model related classes
//#include "GaudiExamples/MyTrack.h"
//
// dd4pod's tracker hit is the input collectiopn
#include "dd4pod/TrackerHitCollection.h"
// edm4hep's tracker hit is the input collectiopn
#include "edm4hep/MCParticle.h"
#include "edm4hep/SimTrackerHitCollection.h"
// eicd's RawTrackerHit is the output
#include "eicd/RawTrackerHitCollection.h"
......@@ -29,7 +28,7 @@ namespace Jug::Digi {
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10.};
Gaudi::Property<double> m_threshold{this, "threshold", 0. * Gaudi::Units::keV};
Rndm::Numbers m_gaussDist;
DataHandle<dd4pod::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
DataHandle<edm4hep::SimTrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};
......@@ -57,7 +56,7 @@ namespace Jug::Digi {
StatusCode execute() override
{
// input collection
const dd4pod::TrackerHitCollection* simhits = m_inputHitCollection.get();
auto simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
// eic::RawTrackerHitCollection* rawHitCollection = new eic::RawTrackerHitCollection();
......@@ -65,35 +64,35 @@ namespace Jug::Digi {
int ID = 0;
for (const auto& ahit : *simhits) {
if (msgLevel(MSG::DEBUG)) {
debug() << "--------------------" << ahit.cellID() << endmsg;
debug() << "Hit in cellID = " << ahit.cellID() << endmsg;
debug() << " position = (" << ahit.position().x << "," << ahit.position().y <<","<< ahit.position().z << ")" << endmsg;
debug() << " xy_radius = " << std::hypot(ahit.position().x , ahit.position().y ) << endmsg;
debug() << " momentum = (" << ahit.momentum().x << "," << ahit.momentum().y <<","<< ahit.momentum().z << ")" << endmsg;
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.energyDeposit() * Gaudi::Units::keV < m_threshold) {
if (ahit.getEDep() * Gaudi::Units::keV < m_threshold) {
if (msgLevel(MSG::DEBUG)) {