Commit 654bdfb4 authored by Whitney Armstrong's avatar Whitney Armstrong

Added new source tracking hit source linker

- TrackingHitsSourceLinker takes multiple TrackerHitCollections
- Creates the source link for all hits (= tracker hit + surface object )
parent 7c19e51b
...@@ -16,110 +16,65 @@ ...@@ -16,110 +16,65 @@
// eicd's RawTrackerHit is the output // eicd's RawTrackerHit is the output
#include "eicd/RawTrackerHitCollection.h" #include "eicd/RawTrackerHitCollection.h"
namespace Jug { namespace Jug::Digi {
namespace Digi {
/** Ultra-fast silicon detector digitization.
*
*
*/
class UFSDTrackerDigi : public GaudiAlgorithm {
public:
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10};
Rndm::Numbers m_gaussDist;
DataHandle<dd4pod::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
public: /** Ultra-fast silicon detector digitization.
*
*
*/
class UFSDTrackerDigi : public GaudiAlgorithm {
public:
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10};
Rndm::Numbers m_gaussDist;
DataHandle<dd4pod::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<eic::RawTrackerHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
this};
public:
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm; // ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
UFSDTrackerDigi(const std::string& name, ISvcLocator* svcLoc) UFSDTrackerDigi(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
: GaudiAlgorithm(name, svcLoc) { {
declareProperty("inputHitCollection", m_inputHitCollection,""); declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, ""); declareProperty("outputHitCollection", m_outputHitCollection, "");
} }
StatusCode initialize() override { StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) if (GaudiAlgorithm::initialize().isFailure())
return StatusCode::FAILURE; return StatusCode::FAILURE;
IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true); IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
StatusCode sc = m_gaussDist.initialize( randSvc, Rndm::Gauss(0.0, m_timeResolution.value())); StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(0.0, m_timeResolution.value()));
if (!sc.isSuccess()) { if (!sc.isSuccess()) {
return StatusCode::FAILURE; return StatusCode::FAILURE;
} }
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
StatusCode execute() override { StatusCode execute() override
{
// input collection // input collection
const dd4pod::TrackerHitCollection* simhits = m_inputHitCollection.get(); const dd4pod::TrackerHitCollection* simhits = m_inputHitCollection.get();
// Create output collections // Create output collections
auto rawhits = m_outputHitCollection.createAndPut(); auto rawhits = m_outputHitCollection.createAndPut();
//eic::RawTrackerHitCollection* rawHitCollection = new eic::RawTrackerHitCollection(); // eic::RawTrackerHitCollection* rawHitCollection = new eic::RawTrackerHitCollection();
std::map<long long,int> cell_hit_map; std::map<long long, int> cell_hit_map;
for(const auto& ahit : *simhits) { for (const auto& ahit : *simhits) {
//std::cout << ahit << "\n"; // std::cout << ahit << "\n";
if(cell_hit_map.count(ahit.cellID()) == 0) { if (cell_hit_map.count(ahit.cellID()) == 0) {
cell_hit_map[ahit.cellID()] = rawhits->size(); cell_hit_map[ahit.cellID()] = rawhits->size();
eic::RawTrackerHit rawhit((long long)ahit.cellID(), eic::RawTrackerHit rawhit((long long)ahit.cellID(),
ahit.truth().time*1e6+m_gaussDist()*1e3, // ns->fs ahit.truth().time * 1e6 + m_gaussDist() * 1e3, // ns->fs
int(ahit.energyDeposit() * 1e6)); int(ahit.energyDeposit() * 1e6));
rawhits->push_back(rawhit); rawhits->push_back(rawhit);
} else { } else {
auto hit = (*rawhits)[cell_hit_map[ahit.cellID()]]; auto hit = (*rawhits)[cell_hit_map[ahit.cellID()]];
hit.time(ahit.truth().time*1e6+m_gaussDist()*1e3); hit.time(ahit.truth().time * 1e6 + m_gaussDist() * 1e3);
auto ch = hit.charge(); auto ch = hit.charge();
hit.charge( ch + int(ahit.energyDeposit() * 1e6)); hit.charge(ch + int(ahit.energyDeposit() * 1e6));
} }
} }
return StatusCode::SUCCESS; return StatusCode::SUCCESS;
} }
}; };
DECLARE_COMPONENT(UFSDTrackerDigi) DECLARE_COMPONENT(UFSDTrackerDigi)
// class DataProducerProp : public GaudiAlgorithm { } // namespace Jug::Digi
// public:
// // ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
// DataProducerProp(const std::string& name, ISvcLocator* svcLoc)
// : GaudiAlgorithm(name, svcLoc) {}
// StatusCode initialize() override {
// f_counter = m_starting_value.value();
// return StatusCode::SUCCESS;
// }
// StatusCode execute() override {
// m_vec.put(ThreeVecEx{f_counter, f_counter + 1, f_counter + 2});
// info() << "executing DataProducer (prop): " << f_counter << " ..."
// << endmsg;
// f_counter++;
// return StatusCode::SUCCESS;
// }
// private:
// Gaudi::Property<int> m_starting_value{this, "StartingValue", 0};
// AnyDataHandle<ThreeVecEx> m_vec{"/Event/UnknownVec2",
// Gaudi::DataHandle::Writer, this};
// int f_counter = 0;
//};
// DECLARE_COMPONENT(DataProducerProp)
// class DataConsumerProp : public GaudiAlgorithm {
// public:
// // using GaudiAlgorithm::GaudiAlgorithm;
// DataConsumerProp(const std::string& name, ISvcLocator* svcLoc)
// : GaudiAlgorithm(name, svcLoc) {}
// StatusCode execute() override {
// info() << "executing DataConsumer (prop): {" << *m_vec.get() << "}"
// << endmsg;
// return StatusCode::SUCCESS;
// }
// private:
// AnyDataHandle<ThreeVecEx> m_vec{"/Event/UnknownVec2",
// Gaudi::DataHandle::Reader, this};
//};
// DECLARE_COMPONENT(DataConsumerProp)
} // namespace Examples
} // namespace Gaudi
...@@ -30,6 +30,7 @@ gaudi_depends_on_subdirs(GaudiAlg GaudiKernel) ...@@ -30,6 +30,7 @@ gaudi_depends_on_subdirs(GaudiAlg GaudiKernel)
gaudi_add_module(JugRecoPlugins gaudi_add_module(JugRecoPlugins
src/components/TrackerHitReconstruction.cpp src/components/TrackerHitReconstruction.cpp
src/components/TrackerSourceLinker.cpp src/components/TrackerSourceLinker.cpp
src/components/TrackingHitsSourceLinker.cpp
src/components/TrackFindingAlgorithm.cpp src/components/TrackFindingAlgorithm.cpp
src/components/TestACTSLogger.cpp src/components/TestACTSLogger.cpp
src/components/TrackParamTruthInit.cpp src/components/TrackParamTruthInit.cpp
......
...@@ -79,15 +79,14 @@ namespace Jug::Reco { ...@@ -79,15 +79,14 @@ namespace Jug::Reco {
using Acts::UnitConstants::ns; using Acts::UnitConstants::ns;
double p_cluster = c.energy()*GeV; double p_cluster = c.energy()*GeV;
if( p_cluster < 1.0) { if( p_cluster/GeV < 0.1) {
debug() << " skipping cluster with energy " << p_cluster/GeV << " GeV" << endmsg; debug() << " skipping cluster with energy " << p_cluster/GeV << " GeV" << endmsg;
continue; continue;
} }
for (const auto& t : *vtx_hits) {
for(const auto& t : *vtx_hits) { double len = std::hypot(t.x(), t.y(), t.z());
double len = std::hypot( t.x() , t.y() , t.z() );
// build some track cov matrix // build some track cov matrix
Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero(); Acts::BoundSymMatrix cov = Acts::BoundSymMatrix::Zero();
...@@ -95,13 +94,16 @@ namespace Jug::Reco { ...@@ -95,13 +94,16 @@ namespace Jug::Reco {
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm; cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = 1.0 * mm*1.0 * mm;
cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0; cov(Acts::eBoundPhi, Acts::eBoundPhi) = M_PI / 180.0;
cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0; cov(Acts::eBoundTheta, Acts::eBoundTheta) = M_PI / 180.0;
cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (0.01*0.01*p_cluster*p_cluster); cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 1.0 / (0.9*0.9*p_cluster*p_cluster);
cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns; cov(Acts::eBoundTime, Acts::eBoundTime) = Acts::UnitConstants::ns;
// add all charges to the track candidate... // add all charges to the track candidate...
init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0), init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
Acts::Vector3D(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len), p_cluster, -1, Acts::Vector3D(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len), p_cluster, -1,
std::make_optional(cov)); std::make_optional(cov));
init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
Acts::Vector3D(t.x() * p_cluster / len, t.y() * p_cluster / len, t.z() * p_cluster / len), p_cluster, 1,
std::make_optional(cov));
} }
//init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0), //init_trk_params->emplace_back(Acts::Vector4D(0 * mm, 0 * mm, 0 * mm, 0),
// Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1, // Acts::Vector3D(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
......
// This file is part of the Acts project.
#include "JugReco/GeometryContainers.hpp"
// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.h"
#include "GaudiKernel/Property.h"
#include "Gaudi/Property.h"
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "DDRec/CellIDPositionConverter.h"
#include "DDRec/SurfaceManager.h"
#include "DDRec/Surface.h"
#include "DD4hep/Volumes.h"
#include "DD4hep/DD4hepUnits.h"
#include "Acts/Utilities/Units.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
#include "JugReco/SourceLinks.h"
#include "eicd/TrackerHitCollection.h"
namespace Jug::Reco {
class TrackingHitsSourceLinker : public GaudiAlgorithm {
public:
using TrkHits = eic::TrackerHitCollection;
Gaudi::Property<std::vector<std::string>> m_inputs{this, "inputTrackerCollections", {}, "List of intput tracker hit collections"};
std::vector<DataHandle<TrkHits>*> m_inputHandles;
// Gaudi::Property<std::vector<std::string>> m_inputCollections{this, "inputCollections", {}, "Input Tracker hit
// collections"}; DataHandle<eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection",
// Gaudi::DataHandle::Reader, this};
DataHandle<SourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
/// Pointer to the geometry service
SmartIF<IGeoSvc> m_geoSvc;
/// Lookup container for hit surfaces that generate smeared hits
std::unordered_map<uint64_t, const Acts::Surface*> m_surfaces;
public:
TrackingHitsSourceLinker(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
{
// declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputSourceLinks", m_outputSourceLinks, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure())
return StatusCode::FAILURE;
for (const auto& col : m_inputs) {
m_inputHandles.push_back(new DataHandle<TrkHits>(col, Gaudi::DataHandle::Reader, this));
// declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("dummy_in_" + col, *(m_inputHandles.back()),"");
}
m_geoSvc = service("GeoSvc");
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;
}
debug() << "visiting all the surfaces " << endmsg;
m_geoSvc->trackingGeometry()->visitSurfaces([this](const Acts::Surface* surface) {
// for now we just require a valid surface
if (not surface) {
return;
}
auto det_element = dynamic_cast<const Acts::DD4hepDetectorElement*>(surface->associatedDetectorElement());
if (!det_element) {
debug() << "invalid det_element!!! " << endmsg;
return;
}
this->m_surfaces.insert_or_assign(det_element->identifier(), surface);
});
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// Create output collections
auto source_links = m_outputSourceLinks.createAndPut();
for (const auto& handle : m_inputHandles) {
// input collection
const eic::TrackerHitCollection* hits = handle->get();
debug() << (*hits).size() << " hits " << endmsg;
for (const auto& ahit : *hits) {
// setup local covariance
Acts::BoundMatrix cov = Acts::BoundMatrix::Zero();
cov(Acts::eBoundLoc0, Acts::eBoundLoc0) =
ahit.covsym_xx() * Acts::UnitConstants::mm * ahit.covsym_xx() * Acts::UnitConstants::mm;
cov(Acts::eBoundLoc1, Acts::eBoundLoc1) =
ahit.covsym_yy() * Acts::UnitConstants::mm * ahit.covsym_yy() * Acts::UnitConstants::mm;
auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
auto vol_id = vol_ctx->identifier;
// debug() << " hit : \n" << ahit << endmsg;
// debug() << "cell ID : " << ahit.cellID() << endmsg;
// debug() << " position : (" << ahit.position(0) << ", " << ahit.position(1) << ", "<< ahit.position(2) <<
// ") " << endmsg;
debug() << " vol_id : " << vol_id << endmsg;
debug() << " placment pos : " << vol_ctx->volumePlacement().position() << endmsg;
const auto is = m_surfaces.find(vol_id);
if (is == m_surfaces.end()) {
debug() << " vol_id (" << vol_id << ") not found in m_surfaces." << endmsg;
continue;
}
const Acts::Surface* surface = is->second;
debug() << " surface center : " << surface->center(Acts::GeometryContext()) << endmsg;
// transform global position into local coordinates
Acts::Vector2D pos(0, 0);
// geometry context contains nothing here
pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.x(), ahit.y(), ahit.z()}, {0, 0, 0})
.value(); //, pos);
//// smear truth to create local measurement
Acts::BoundVector loc = Acts::BoundVector::Zero();
loc[Acts::eBoundLoc0] = pos[0]; //+ m_cfg.sigmaLoc0 * stdNormal(rng);
loc[Acts::eBoundLoc0] = pos[1]; //+ m_cfg.sigmaLoc1 * stdNormal(rng);
debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
// create source link at the end of the container
// auto it = source_links->emplace_hint(source_links->end(), *surface, hit, 2, loc, cov);
auto it = source_links->emplace_hint(source_links->end(), *surface, 2, loc, cov);
// ensure hits and links share the same order to prevent ugly surprises
if (std::next(it) != source_links->end()) {
error() << "The hit ordering broke. Run for your life." << endmsg;
return StatusCode::FAILURE;
}
}
}
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(TrackingHitsSourceLinker)
} // namespace Jug::Reco
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment