198 larprop = lar::providerFrom<detinfo::LArPropertiesService>();
201 std::vector<std::map<int, double>> vtimemap(3);
202 std::vector<std::map<int, art::Ptr<recob::Hit>>> vhitmap(3);
205 for (
size_t ihit = 0; ihit < fHits.size(); ++ihit) {
207 unsigned int w = hitWireID.
Wire;
208 unsigned int pl = hitWireID.
Plane;
209 double time = fHits[ihit]->PeakTime();
210 time -= detProp.GetXTicksOffset(
212 vtimemap[pl][
w] = time;
213 vhitmap[pl][
w] = fHits[ihit];
220 unsigned maxnumhits0 = 0;
221 unsigned maxnumhits1 = 0;
223 std::vector<double> tmin(vtimemap.size());
224 std::vector<double> tmax(vtimemap.size());
225 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
230 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
231 for (
auto itime = vtimemap[iclu].
begin(); itime != vtimemap[iclu].end(); ++itime) {
232 if (itime->second > tmax[iclu]) { tmax[iclu] = itime->second; }
233 if (itime->second < tmin[iclu]) { tmin[iclu] = itime->second; }
235 if (vtimemap[iclu].
size() > maxnumhits0) {
238 maxnumhits1 = maxnumhits0;
241 maxnumhits0 = vtimemap[iclu].size();
243 else if (vtimemap[iclu].
size() > maxnumhits1) {
245 maxnumhits1 = vtimemap[iclu].size();
252 for (
int iclu = 0; iclu < (int)vtimemap.size(); ++iclu) {
253 if (iclu != iclu1 && iclu != iclu2) iclu3 = iclu;
256 if (iclu1 != -1 && iclu2 != -1) {
258 auto ihit = vhitmap[iclu1].begin();
259 auto itime = vtimemap[iclu1].begin();
260 while (itime != vtimemap[iclu1].
end()) {
261 if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) -
ftmatch ||
262 itime->second > std::min(tmax[iclu1], tmax[iclu2]) +
ftmatch) {
263 vtimemap[iclu1].erase(itime++);
264 vhitmap[iclu1].erase(ihit++);
272 ihit = vhitmap[iclu2].begin();
273 itime = vtimemap[iclu2].begin();
274 while (itime != vtimemap[iclu2].
end()) {
275 if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) -
ftmatch ||
276 itime->second > std::min(tmax[iclu1], tmax[iclu2]) +
ftmatch) {
277 vtimemap[iclu2].erase(itime++);
278 vhitmap[iclu2].erase(ihit++);
287 if (!vtimemap[iclu1].
size()) {
288 if (iclu3 != -1) {
std::swap(iclu3, iclu1); }
290 if (!vtimemap[iclu2].
size()) {
296 if ((!vtimemap[iclu1].
size()) || (!vtimemap[iclu2].
size()))
return;
299 auto times1 = vtimemap[iclu1].begin();
300 auto timee1 = vtimemap[iclu1].end();
302 auto times2 = vtimemap[iclu2].begin();
303 auto timee2 = vtimemap[iclu2].end();
306 double ts1 = times1->second;
307 double te1 = timee1->second;
308 double ts2 = times2->second;
309 double te2 = timee2->second;
317 double Efield_drift = detProp.Efield(0);
318 double Temperature = detProp.Temperature();
320 double driftvelocity = detProp.DriftVelocity(
321 Efield_drift, Temperature);
322 double timepitch = driftvelocity * timetick;
328 std::vector<double> vtracklength;
329 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
331 double tracklength = 0;
332 if (vtimemap[iclu].
size() == 1) { tracklength = wire_pitch; }
333 else if (!vtimemap[iclu].
empty()) {
335 wend = vtimemap[iclu].cend();
336 double t0 = iw->first, w0 = iw->second;
337 while (++iw != wend) {
338 tracklength += std::hypot((iw->first - w0) * wire_pitch, (iw->second - t0) * timepitch);
343 vtracklength.push_back(tracklength);
346 std::map<int, int> maxhitsMatch;
348 auto ihit1 = vhitmap[iclu1].begin();
349 for (
auto itime1 = vtimemap[iclu1].
begin(); itime1 != vtimemap[iclu1].end();
351 std::vector<art::Ptr<recob::Hit>> sp_hits;
352 sp_hits.push_back(ihit1->second);
353 double hitcoord[3] = {0., 0., 0.};
355 if (vtimemap[iclu1].
size() == 1) { length1 = wire_pitch; }
357 for (
auto iw1 = vtimemap[iclu1].
begin(); iw1 != itime1; ++iw1) {
360 length1 += std::hypot((iw1->first - iw2->first) * wire_pitch,
361 (iw1->second - iw2->second) * timepitch);
364 double difference = 1e10;
365 auto matchedtime = vtimemap[iclu2].end();
366 auto matchedhit = vhitmap[iclu2].end();
368 auto ihit2 = vhitmap[iclu2].begin();
369 for (
auto itime2 = vtimemap[iclu2].
begin(); itime2 != vtimemap[iclu2].end();
371 if (maxhitsMatch[itime2->first])
continue;
373 if (vtimemap[iclu2].
size() == 1) { length2 = wire_pitch; }
375 for (
auto iw1 = vtimemap[iclu2].
begin(); iw1 != itime2; ++iw1) {
378 length2 += std::hypot((iw1->first - iw2->first) * wire_pitch,
379 (iw1->second - iw2->second) * timepitch);
382 if (rev) length2 = vtracklength[iclu2] - length2;
383 length2 = vtracklength[iclu1] / vtracklength[iclu2] * length2;
385 if (timematch &&
std::abs(length2 - length1) < difference) {
386 difference =
std::abs(length2 - length1);
387 matchedtime = itime2;
392 hitcoord[0] = detProp.ConvertTicksToX(
393 matchedtime->second + detProp.GetXTicksOffset((matchedhit->second)->WireID().Plane,
394 (matchedhit->second)->
WireID().TPC,
395 (matchedhit->second)->
WireID().Cryostat),
396 (matchedhit->second)->WireID().Plane,
397 (matchedhit->second)->
WireID().TPC,
398 (matchedhit->second)->
WireID().Cryostat);
409 hitcoord[1] = intersection->y;
410 hitcoord[2] = intersection->z;
413 if (hitcoord[1] > -1e9 && hitcoord[2] > -1e9) {
414 maxhitsMatch[matchedtime->first] = 1;
415 sp_hits.push_back(matchedhit->second);
418 if (sp_hits.size() > 1) {
419 trajPos.push_back(TVector3(hitcoord));
code to link reconstructed objects back to the MC truth information
double ftmatch
tolerance for time matching (in ticks)
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)
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
geo::WireReadoutGeom const * wireReadoutGeom
const detinfo::LArProperties * larprop
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.
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
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
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
std::vector< std::vector< art::Ptr< recob::Hit > > > trajHit