Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
#include "JugTrack/GeometryContainers.hpp"
// Gaudi
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiAlg/Transformer.h"
#include "GaudiAlg/GaudiTool.h"
#include "GaudiKernel/RndmGenerators.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/Surfaces/Surface.hpp"
#include "Acts/Definitions/Units.hpp"
#include "Acts/Definitions/Common.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
#include "JugTrack/Index.hpp"
#include "JugTrack/IndexSourceLink.hpp"
#include "JugTrack/Measurement.hpp"
#include "JugTrack/ProtoTrack.hpp"
#include "eicd/TrackerHitCollection.h"
namespace Jug::Reco {
/** Single Track source Linker and proto tracks.
* This algorithm assumes only single track events.
*
*/
class SingleTrackSourceLinker : public GaudiAlgorithm {
public:
DataHandle<eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<IndexSourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
DataHandle<MeasurementContainer> m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this};
DataHandle<ProtoTrackContainer> m_outputProtoTracks{"outputProtoTracks", 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:
SingleTrackSourceLinker(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc) {
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputSourceLinks", m_outputSourceLinks, "");
declareProperty("outputMeasurements", m_outputMeasurements, "");
declareProperty("outputProtoTracks", m_outputProtoTracks, "");
}
StatusCode initialize() override {
if (GaudiAlgorithm::initialize().isFailure())
return StatusCode::FAILURE;
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;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override {
// input collection
const eic::TrackerHitCollection* hits = m_inputHitCollection.get();
// Create output collections
auto sourceLinks = m_outputSourceLinks.createAndPut();
auto measurements = m_outputMeasurements.createAndPut();
auto protoTracks = m_outputProtoTracks.createAndPut();
// IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
// IndexMultimap<Index> hitSimHitsMap;
sourceLinks->reserve(hits->size());
measurements->reserve(hits->size());
// assume single track --> one ProtoTrack
ProtoTrack track;
track.reserve((*hits).size());
debug() << (*hits).size() << " hits " << endmsg;
int ihit = 0;
for(const auto& ahit : *hits) {
track.emplace_back(ihit);
auto vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.cellID());
auto vol_id = vol_ctx->identifier;
const auto is = m_geoSvc->surfaceMap().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;
// NOTE
// Here is where the important hit and tracking geometry is connected.
// A "Measurement" is constructed to for each hit which makes the connection to
// the tracking surface and covariance matrix
auto volman = m_geoSvc->detector()->volumeManager();
auto detelem = volman.lookupDetElement(vol_id);
auto local_pos = detelem.nominal().worldToLocal({ahit.x(), ahit.y(), ahit.z()});
auto acts_pos = surface->globalToLocal(Acts::GeometryContext(), {ahit.x(), ahit.y(), ahit.z()}, {0, 0, 0}).value();
debug() << "===== Debugging hit =====" << endmsg;
debug() << "DD4hep global pos (" << ahit.x() << "," << ahit.y() << "," << ahit.z() << ")" << endmsg;
debug() << "DD4hep local pos (" << local_pos.x() << "," << local_pos.y() << "," << local_pos.z() << ")" << endmsg;
debug() << "ACTS local position : (" << acts_pos[0] << "," << acts_pos[1] << ")"<< endmsg;
debug() << "ACTS surface center : " << surface->center(Acts::GeometryContext()).transpose() << endmsg;
debug() << "DD4hep DetElement center : " << detelem.nominal().localToWorld(detelem.placement().position())/dd4hep::mm << endmsg;
// construct the vector of measured parameters (2d position in this case)
Acts::Vector2 pos(acts_pos.x(), acts_pos.y());
// construct the covariance matrix
Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
cov(0,0) = ahit.covsym_xx()*Acts::UnitConstants::mm*ahit.covsym_xx()*Acts::UnitConstants::mm;
cov(1,1) = ahit.covsym_yy()*Acts::UnitConstants::mm*ahit.covsym_yy()*Acts::UnitConstants::mm;
// Above we only consider the two position coordinates the comment below shows how to add time
// which we will probably want to try later.
//
//Acts::SymMatrix3 cov;
//cov << 0.05, 0., 0., 0., 0.05, 0., 0., 0., 900. * Acts::UnitConstants::ps * Acts::UnitConstants::ps;
//Acts::Vector3 par(localX, localY, simHit.time());
Index hitIdx = ihit;
IndexSourceLink sourceLink(surface->geometryId(), ihit);
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
auto meas = Acts::makeMeasurement(sourceLink, pos, cov,
Acts::eBoundLoc0, Acts::eBoundLoc1);//, Acts::eBoundTime);
// TODO: check that local to global is the same in dd4hep and acts.
//
//debug() << "ACTS surface center : " << surface->center(Acts::GeometryContext()) << endmsg;
// transform global position into local coordinates
// geometry context contains nothing here
//Acts::Vector2 loc = Acts::Vector2::Zero();
//loc[Acts::eBoundLoc0] = pos[0] ;//+ m_cfg.sigmaLoc0 * stdNormal(rng);
//loc[Acts::eBoundLoc1] = pos[1] ;//+ m_cfg.sigmaLoc1 * stdNormal(rng);
////debug() << "loc : (" << loc[0] << ", " << loc[1] << ")" << endmsg;
////local position
////auto loc = {ahit.x(), ahit.y(), ahit.z()} - vol_ctx->volumePlacement().position()
////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;
//// the measurement container is unordered and the index under which the
//// measurement will be stored is known before adding it.
//auto meas = Acts::makeMeasurement(sourceLink, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
// add to output containers. since the input is already geometry-order,
// new elements in geometry containers can just be appended at the end.
sourceLinks->emplace_hint(sourceLinks->end(), std::move(sourceLink));
measurements->emplace_back(std::move(meas));
ihit++;
}
// add proto track to the output collection
protoTracks->emplace_back(std::move(track));
return StatusCode::SUCCESS;
}
};
DECLARE_COMPONENT(SingleTrackSourceLinker)
} // namespace Jug::reco