196 larprop = lar::providerFrom<detinfo::LArPropertiesService>();
199 std::vector<std::map<int, double>> vtimemap(3);
200 std::vector<std::map<int, art::Ptr<recob::Hit>>> vhitmap(3);
203 for (
size_t ihit = 0; ihit < fHits.size(); ++ihit) {
205 unsigned int w = hitWireID.
Wire;
206 unsigned int pl = hitWireID.
Plane;
207 double time = fHits[ihit]->PeakTime();
208 time -= detProp.GetXTicksOffset(
210 vtimemap[pl][
w] = time;
211 vhitmap[pl][
w] = fHits[ihit];
218 unsigned maxnumhits0 = 0;
219 unsigned maxnumhits1 = 0;
221 std::vector<double> tmin(vtimemap.size());
222 std::vector<double> tmax(vtimemap.size());
223 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
228 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
229 for (
auto itime = vtimemap[iclu].
begin(); itime != vtimemap[iclu].end(); ++itime) {
230 if (itime->second > tmax[iclu]) { tmax[iclu] = itime->second; }
231 if (itime->second < tmin[iclu]) { tmin[iclu] = itime->second; }
233 if (vtimemap[iclu].
size() > maxnumhits0) {
236 maxnumhits1 = maxnumhits0;
239 maxnumhits0 = vtimemap[iclu].size();
241 else if (vtimemap[iclu].
size() > maxnumhits1) {
243 maxnumhits1 = vtimemap[iclu].size();
250 for (
int iclu = 0; iclu < (int)vtimemap.size(); ++iclu) {
251 if (iclu != iclu1 && iclu != iclu2) iclu3 = iclu;
254 if (iclu1 != -1 && iclu2 != -1) {
256 auto ihit = vhitmap[iclu1].begin();
257 auto itime = vtimemap[iclu1].begin();
258 while (itime != vtimemap[iclu1].
end()) {
259 if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) -
ftmatch ||
260 itime->second > std::min(tmax[iclu1], tmax[iclu2]) +
ftmatch) {
261 vtimemap[iclu1].erase(itime++);
262 vhitmap[iclu1].erase(ihit++);
270 ihit = vhitmap[iclu2].begin();
271 itime = vtimemap[iclu2].begin();
272 while (itime != vtimemap[iclu2].
end()) {
273 if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) -
ftmatch ||
274 itime->second > std::min(tmax[iclu1], tmax[iclu2]) +
ftmatch) {
275 vtimemap[iclu2].erase(itime++);
276 vhitmap[iclu2].erase(ihit++);
285 if (!vtimemap[iclu1].
size()) {
286 if (iclu3 != -1) {
std::swap(iclu3, iclu1); }
288 if (!vtimemap[iclu2].
size()) {
294 if ((!vtimemap[iclu1].
size()) || (!vtimemap[iclu2].
size()))
return;
297 auto times1 = vtimemap[iclu1].begin();
298 auto timee1 = vtimemap[iclu1].end();
300 auto times2 = vtimemap[iclu2].begin();
301 auto timee2 = vtimemap[iclu2].end();
304 double ts1 = times1->second;
305 double te1 = timee1->second;
306 double ts2 = times2->second;
307 double te2 = timee2->second;
315 double Efield_drift = detProp.Efield(0);
316 double Temperature = detProp.Temperature();
318 double driftvelocity = detProp.DriftVelocity(
319 Efield_drift, Temperature);
320 double timepitch = driftvelocity * timetick;
326 std::vector<double> vtracklength;
327 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
329 double tracklength = 0;
330 if (vtimemap[iclu].
size() == 1) { tracklength = wire_pitch; }
331 else if (!vtimemap[iclu].
empty()) {
333 wend = vtimemap[iclu].cend();
334 double t0 = iw->first, w0 = iw->second;
335 while (++iw != wend) {
336 tracklength += std::hypot((iw->first - w0) * wire_pitch, (iw->second - t0) * timepitch);
341 vtracklength.push_back(tracklength);
344 std::map<int, int> maxhitsMatch;
346 auto ihit1 = vhitmap[iclu1].begin();
347 for (
auto itime1 = vtimemap[iclu1].
begin(); itime1 != vtimemap[iclu1].end();
349 std::vector<art::Ptr<recob::Hit>> sp_hits;
350 sp_hits.push_back(ihit1->second);
351 double hitcoord[3] = {0., 0., 0.};
353 if (vtimemap[iclu1].
size() == 1) { length1 = wire_pitch; }
355 for (
auto iw1 = vtimemap[iclu1].
begin(); iw1 != itime1; ++iw1) {
358 length1 += std::hypot((iw1->first - iw2->first) * wire_pitch,
359 (iw1->second - iw2->second) * timepitch);
362 double difference = 1e10;
363 auto matchedtime = vtimemap[iclu2].end();
364 auto matchedhit = vhitmap[iclu2].end();
366 auto ihit2 = vhitmap[iclu2].begin();
367 for (
auto itime2 = vtimemap[iclu2].
begin(); itime2 != vtimemap[iclu2].end();
369 if (maxhitsMatch[itime2->first])
continue;
371 if (vtimemap[iclu2].
size() == 1) { length2 = wire_pitch; }
373 for (
auto iw1 = vtimemap[iclu2].
begin(); iw1 != itime2; ++iw1) {
376 length2 += std::hypot((iw1->first - iw2->first) * wire_pitch,
377 (iw1->second - iw2->second) * timepitch);
380 if (rev) length2 = vtracklength[iclu2] - length2;
381 length2 = vtracklength[iclu1] / vtracklength[iclu2] * length2;
383 if (timematch &&
std::abs(length2 - length1) < difference) {
384 difference =
std::abs(length2 - length1);
385 matchedtime = itime2;
390 hitcoord[0] = detProp.ConvertTicksToX(
391 matchedtime->second + detProp.GetXTicksOffset((matchedhit->second)->WireID().Plane,
392 (matchedhit->second)->
WireID().TPC,
393 (matchedhit->second)->
WireID().Cryostat),
394 (matchedhit->second)->WireID().Plane,
395 (matchedhit->second)->
WireID().TPC,
396 (matchedhit->second)->
WireID().Cryostat);
410 hitcoord[1] = tmpWIDI.
y;
411 hitcoord[2] = tmpWIDI.
z;
414 if (hitcoord[1] > -1e9 && hitcoord[2] > -1e9) {
415 maxhitsMatch[matchedtime->first] = 1;
416 sp_hits.push_back(matchedhit->second);
419 if (sp_hits.size() > 1) {
420 trajPos.push_back(TVector3(hitcoord));
code to link reconstructed objects back to the MC truth information
double z
z position of intersection
double ftmatch
tolerance for time matching (in ticks)
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< TVector3 > trajPos
WireID_t Wire
Index of the wire within its plane.
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
PlaneID_t Plane
Index of the plane within its TPC.
double fsmatch
tolerance for distance matching (in cm)
art::ServiceHandle< geo::Geometry const > geom
const detinfo::LArProperties * larprop
double y
y position of intersection
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
second_as<> second
Type of time stored in seconds, in double precision.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
std::vector< std::vector< art::Ptr< recob::Hit > > > trajHit