Skip to content

GitLab

Commit 4fd6165b authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

Eicd update

parent ece25573
Hide whitespace changes
Inline Side-by-side
JugBase/CMakeLists.txt
View file @ 4fd6165b
......@@ -34,7 +34,6 @@ gaudi_add_module(JugBasePlugins
src/components/GeoSvc.cpp
src/components/ParticleSvc.cpp
src/components/InputCopier.cpp
src/components/MC2DummyParticle.cpp
src/components/PodioInput.cpp
src/components/PodioOutput.cpp
src/components/ReadTestConsumer.cxx
......
JugBase/src/components/MC2DummyParticle.cpp deleted 100644 → 0
View file @ ece25573
#include "Gaudi/Algorithm.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Producer.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/RndmGenerators.h"
#include <algorithm>
#include <cmath>
#include "JugBase/DataHandle.h"
// Event Model related classes
#include "edm4hep/MCParticleCollection.h"
#include "eicd/ReconstructedParticleCollection.h"
namespace Jug::Base {
class MC2DummyParticle : public GaudiAlgorithm {
public:
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*/};
MC2DummyParticle(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("inputCollection", m_inputParticles, "MCParticles");
declareProperty("outputCollection", m_outputParticles, "DummyReconstructedParticles");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure())
return StatusCode::FAILURE;
IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(1.0, m_smearing.value()));
if (!sc.isSuccess()) {
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collection
const edm4hep::MCParticleCollection* parts = m_inputParticles.get();
// output collection
auto& out_parts = *(m_outputParticles.createAndPut());
for (const auto& p : *parts) {
if (p.getGeneratorStatus() > 1) {
if (msgLevel(MSG::DEBUG)) {
debug() << "ignoring particle with generatorStatus = " << p.getGeneratorStatus() << endmsg;
}
continue;
}
// 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 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 auto vx = p.getVertex().x;
const auto vy = p.getVertex().y;
const auto vz = p.getVertex().z;
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.p({px, py, pz});
rec_part.v({vx, vy, vz});
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.theta(p_theta);
rec_part.phi(p_phi);
rec_part.momentum(momentum);
rec_part.energy(energy);
rec_part.mass(p.getMass());
}
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(MC2DummyParticle)
} // namespace Jug::Base
JugDigi/CMakeLists.txt
View file @ 4fd6165b
......@@ -8,7 +8,6 @@ gaudi_add_module(JugDigiPlugins
src/components/CalorimeterHitDigi.cpp
src/components/CalorimeterBirksCorr.cpp
src/components/PhotoMultiplierDigi.cpp
src/components/UFSDTrackerDigi.cpp
src/components/SiliconTrackerDigi.cpp
src/components/SimTrackerHitsCollector.cpp
LINK
......
JugDigi/src/components/CalorimeterHitDigi.cpp
View file @ 4fd6165b
......@@ -77,7 +77,7 @@ namespace Jug::Digi {
DataHandle<edm4hep::SimCalorimeterHitCollection> m_inputHitCollection{
"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawCalorimeterHitCollection> m_outputHitCollection{
DataHandle<eicd::RawCalorimeterHitCollection> m_outputHitCollection{
"outputHitCollection", Gaudi::DataHandle::Writer, this};
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
......@@ -185,7 +185,7 @@ namespace Jug::Digi {
time = c.getTime();
const long long tdc = std::llround((time + m_normDist() * tRes) * stepTDC);
eic::RawCalorimeterHit rawhit(
eicd::RawCalorimeterHit rawhit(
ahit.getCellID(),
(adc > m_capADC.value() ? m_capADC.value() : adc),
tdc
......@@ -238,7 +238,7 @@ namespace Jug::Digi {
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(
eicd::RawCalorimeterHit rawhit(
id,
(adc > m_capADC.value() ? m_capADC.value() : adc),
tdc
......
JugDigi/src/components/PhotoMultiplierDigi.cpp
View file @ 4fd6165b
......@@ -38,7 +38,7 @@ class PhotoMultiplierDigi : public GaudiAlgorithm
public:
DataHandle<edm4hep::SimTrackerHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawPMTHitCollection>
DataHandle<eicd::RawPMTHitCollection>
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}}};
......@@ -121,8 +121,7 @@ public:
// build hit
for (auto &it : hit_groups) {
for (auto &data : it.second) {
eic::RawPMTHit hit{
0 /* Old ID TBDeleted */,
eicd::RawPMTHit hit{
it.first,
static_cast<uint32_t>(data.signal),
static_cast<uint32_t>(data.time/(m_timeStep/ns))};
......
JugDigi/src/components/SiliconTrackerDigi.cpp
View file @ 4fd6165b
#include <algorithm>
#include <cmath>
#include "GaudiAlg/Transformer.h"
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/RndmGenerators.h"
#include "Gaudi/Property.h"
#include "JugBase/DataHandle.h"
......@@ -18,82 +18,80 @@
namespace Jug::Digi {
/** Silicon detector digitization.
*
* \ingroup digi
*/
class SiliconTrackerDigi : public GaudiAlgorithm {
public:
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,
/** Silicon detector digitization.
*
* \ingroup digi
*/
class SiliconTrackerDigi : public GaudiAlgorithm {
public:
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<eicd::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};
DataHandle<eic::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, "");
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;
}
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
auto simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
// eicd::RawTrackerHitCollection* rawHitCollection = new eicd::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;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collection
auto simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
// eic::RawTrackerHitCollection* rawHitCollection = new eic::RawTrackerHitCollection();
std::map<long long, int> cell_hit_map;
for (const auto& ahit : *simhits) {
if (ahit.getEDep() * Gaudi::Units::keV < m_threshold) {
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;
}
debug() << " edep = " << ahit.getEDep() << " (below threshold of " << m_threshold / Gaudi::Units::keV
<< " keV)" << endmsg;
}
if (cell_hit_map.count(ahit.getCellID()) == 0) {
cell_hit_map[ahit.getCellID()] = rawhits->size();
eic::RawTrackerHit rawhit(0 /* deleteme */,
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.getCellID()]];
hit.time(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
auto ch = hit.charge();
hit.charge(ch + std::llround(ahit.getEDep() * 1e6));
continue;
} else {
if (msgLevel(MSG::DEBUG)) {
debug() << " edep = " << ahit.getEDep() << endmsg;
}
}
return StatusCode::SUCCESS;
if (cell_hit_map.count(ahit.getCellID()) == 0) {
cell_hit_map[ahit.getCellID()] = rawhits->size();
eicd::RawTrackerHit rawhit(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.getCellID()]];
hit.timeStamp(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
auto ch = hit.charge();
hit.charge(ch + std::llround(ahit.getEDep() * 1e6));
}
}
};
DECLARE_COMPONENT(SiliconTrackerDigi)
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(SiliconTrackerDigi)
} // namespace Jug::Digi
JugDigi/src/components/UFSDTrackerDigi.cpp deleted 100644 → 0
View file @ ece25573
#include <algorithm>
#include <cmath>
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/RndmGenerators.h"
#include "Gaudi/Property.h"
#include "JugBase/DataHandle.h"
// Event Model related classes
// 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"
namespace Jug::Digi {
/** UFSD detector digitization.
*
* \ingroup digi
*/
class UFSDTrackerDigi : public GaudiAlgorithm {
public:
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10.};
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<eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};
public:
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
UFSDTrackerDigi(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
auto simhits = m_inputHitCollection.get();
// Create output collections
auto rawhits = m_outputHitCollection.createAndPut();
// eic::RawTrackerHitCollection* rawHitCollection = new eic::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;
}
}
// std::cout << ahit << "\n";
if (cell_hit_map.count(ahit.getCellID()) == 0) {
cell_hit_map[ahit.getCellID()] = rawhits->size();
eic::RawTrackerHit rawhit(0 /* TBDeleted */,
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.getCellID()]];
hit.time(ahit.getMCParticle().getTime() * 1e6 + m_gaussDist() * 1e3);
auto ch = hit.charge();
hit.charge(ch + std::llround(ahit.getEDep() * 1e6));
}
}
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(UFSDTrackerDigi)
} // namespace Jug::Digi
JugFast/src/components/ClusterMerger.cpp
View file @ 4fd6165b
......@@ -28,9 +28,9 @@ namespace Jug::Fast {
class ClusterMerger : public GaudiAlgorithm {
public:
// Input
DataHandle<eic::ClusterCollection> m_inputClusters{"InputClusters", Gaudi::DataHandle::Reader, this};
DataHandle<eicd::ClusterCollection> m_inputClusters{"InputClusters", Gaudi::DataHandle::Reader, this};
// Output
DataHandle<eic::ClusterCollection> m_outputClusters{"OutputClusters", Gaudi::DataHandle::Writer, this};
DataHandle<eicd::ClusterCollection> m_outputClusters{"OutputClusters", Gaudi::DataHandle::Writer, this};
public:
ClusterMerger(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
......@@ -84,7 +84,7 @@ public:
float energyError = 0;
float time = 0;
int nhits = 0;
eic::VectorXYZ position;
eicd::VectorXYZ position;
for (const auto& clus : clusters) {
if (msgLevel(MSG::DEBUG)) {
debug() << " --> Adding cluster with energy: " << clus.energy() << endmsg;
......@@ -116,8 +116,8 @@ public:
}
// get a map of mcID --> std::vector<cluster> for clusters that belong together
std::map<eic::Index, std::vector<eic::ConstCluster>> indexedClusterLists(const eic::ClusterCollection& clusters) const {
std::map<eic::Index, std::vector<eic::ConstCluster>> matched = {};
std::map<eicd::Index, std::vector<eicd::ConstCluster>> indexedClusterLists(const eicd::ClusterCollection& clusters) const {
std::map<eicd::Index, std::vector<eicd::ConstCluster>> matched = {};
for (const auto& cluster : clusters) {
if (msgLevel(MSG::VERBOSE)) {
verbose() << " --> Found cluster with mcID " << cluster.mcID() << " and energy "
......
JugFast/src/components/InclusiveKinematicsTruth.cpp
View file @ 4fd6165b
......@@ -21,7 +21,7 @@ class InclusiveKinematicsTruth : public GaudiAlgorithm {
public:
DataHandle<edm4hep::MCParticleCollection> m_inputParticleCollection{"MCParticles", Gaudi::DataHandle::Reader,
this};
DataHandle<eic::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{"InclusiveKinematicsTruth",
DataHandle<eicd::InclusiveKinematicsCollection> m_outputInclusiveKinematicsCollection{"InclusiveKinematicsTruth",
Gaudi::DataHandle::Writer, this};
SmartIF<IParticleSvc> m_pidSvc;
......@@ -139,7 +139,8 @@ public:
kin.nu(q * pi / m_proton);
kin.x(kin.Q2() / (2. * q * pi));
kin.W(sqrt((pi + q).m2()));
kin.scatID(mcscatID);
// kin.scat(????); @TODO: this is now set as a OneToOneRelation to ReconstructedParticle,
// which breaks for this algorithm that takes raw MCParticles
// Debugging output
if (msgLevel(MSG::DEBUG)) {
......
JugFast/src/components/MC2SmearedParticle.cpp
View file @ 4fd6165b
......@@ -17,13 +17,12 @@ namespace Jug::Fast {
class MC2SmearedParticle : public GaudiAlgorithm {
public:
DataHandle<edm4hep::MCParticleCollection> m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ReconstructedParticleCollection> m_outputParticles{"SmearedReconstructedParticles",
Gaudi::DataHandle::Writer, this};
DataHandle<eicd::ReconstructedParticleCollection> m_outputParticles{"SmearedReconstructedParticles",
Gaudi::DataHandle::Writer, this};
Rndm::Numbers m_gaussDist;
Gaudi::Property<double> m_smearing{this, "smearing", 0.01 /* 1 percent*/};
MC2SmearedParticle(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
MC2SmearedParticle(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
declareProperty("inputParticles", m_inputMCParticles, "MCParticles");
declareProperty("outputParticles", m_outputParticles, "SmearedReconstructedParticles");
}
......@@ -54,36 +53,36 @@ public:
// 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 auto& pvec = p.getMomentum();
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();
// momentum smearing
const auto px = p.getMomentum().x * momentum / pgen;
const auto py = p.getMomentum().y * momentum / pgen;
const auto pz = p.getMomentum().z * momentum / pgen;
const auto phi = std::atan2(py, px);
const auto theta = std::atan2(std::hypot(px, py), pz);
// make sure we keep energy consistent
using MomType = decltype(eicd::ReconstructedParticle().momentum().x);
const MomType energy = std::sqrt(p.getEnergy() * p.getEnergy() - pgen * pgen + momentum * momentum);
const MomType px = p.getMomentum().x * momentum / pgen;
const MomType py = p.getMomentum().y * momentum / pgen;
const MomType pz = p.getMomentum().z * momentum / pgen;
const MomType dpx = m_smearing.value() * px;
const MomType dpy = m_smearing.value() * py;
const MomType dpz = m_smearing.value() * pz;
const MomType dE = m_smearing.value() * energy;
// ignore covariance for now
// @TODO: vertex smearing
const auto vx = p.getVertex().x;
const auto vy = p.getVertex().y;
const auto vz = p.getVertex().z;
const MomType vx = p.getVertex().x;
const MomType vy = p.getVertex().y;
const MomType vz = p.getVertex().z;
auto rec_part = out_parts.create();
rec_part.p({px, py, pz});
rec_part.v({vx, vy, vz});
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.theta(theta);
rec_part.phi(phi);
rec_part.momentum(momentum);
rec_part.type(-1); // @TODO: determine type codes
rec_part.energy(energy);
rec_part.momentum({px, py, pz});
rec_part.referencePoint({vx, vy, vz}); // @FIXME: probably not what we want?
rec_part.charge(p.getCharge());
rec_part.mass(p.getMass());
// FIXME disabled mcID link
//rec_part.mcID({p.id(), kMonteCarloSource});
rec_part.goodnessOfPID(1); // Perfect PID
rec_part.covMatrix({dpx, dpy, dpz, dE});
rec_part.PDG(p.getPDG());
}
return StatusCode::SUCCESS;
}
......
JugFast/src/components/MatchClusters.cpp
View file @ 4fd6165b
......@@ -30,13 +30,13 @@ class MatchClusters : public GaudiAlgorithm {
public:
// input data
DataHandle<edm4hep::MCParticleCollection> m_inputMCParticles{"MCParticles", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ReconstructedParticleCollection> m_inputParticles{"ReconstructedChargedParticles",
DataHandle<eicd::ReconstructedParticleCollection> m_inputParticles{"ReconstructedChargedParticles",
Gaudi::DataHandle::Reader, this};
Gaudi::Property<std::vector<std::string>> m_inputClusters{this, "inputClusters", {}, "Clusters to be aggregated"};
std::vector<DataHandle<eic::ClusterCollection>*> m_inputClusterCollections;
std::vector<DataHandle<eicd::ClusterCollection>*> m_inputClusterCollections;
// output data
DataHandle<eic::ReconstructedParticleCollection> m_outputParticles{"ReconstructedParticles",
DataHandle<eicd::ReconstructedParticleCollection> m_outputParticles{"ReconstructedParticles",
Gaudi::DataHandle::Writer, this};
MatchClusters(const std::string& name, ISvcLocator* svcLoc)
......@@ -130,20 +130,20 @@ public:
}
private:
std::vector<DataHandle<eic::ClusterCollection>*> getClusterCollections(const std::vector<std::string>& cols) {
std::vector<DataHandle<eic::ClusterCollection>*> ret;
std::vector<DataHandle<eicd::ClusterCollection>*> getClusterCollections(const std::vector<std::string>& cols) {
std::vector<DataHandle<eicd::ClusterCollection>*> ret;
for (const auto& colname : cols) {
debug() << "initializing cluster collection: " << colname << endmsg;
ret.push_back(new DataHandle<eic::ClusterCollection>{colname, Gaudi::DataHandle::Reader, this});