199 larprop = lar::providerFrom<detinfo::LArPropertiesService>();
202 std::vector<std::map<int, double>> vtimemap(3);
203 std::vector<std::map<int, art::Ptr<recob::Hit>>> vhitmap(3);
206 for (
size_t ihit = 0; ihit < fHits.size(); ++ihit) {
208 unsigned int w = hitWireID.
Wire;
209 unsigned int pl = hitWireID.
Plane;
210 double time = fHits[ihit]->PeakTime();
211 time -= detProp.GetXTicksOffset(
213 vtimemap[pl][
w] = time;
214 vhitmap[pl][
w] = fHits[ihit];
221 unsigned maxnumhits0 = 0;
222 unsigned maxnumhits1 = 0;
224 std::vector<double> tmin(vtimemap.size());
225 std::vector<double> tmax(vtimemap.size());
226 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
231 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
232 for (
auto itime = vtimemap[iclu].
begin(); itime != vtimemap[iclu].end(); ++itime) {
233 if (itime->second > tmax[iclu]) { tmax[iclu] = itime->second; }
234 if (itime->second < tmin[iclu]) { tmin[iclu] = itime->second; }
236 if (vtimemap[iclu].
size() > maxnumhits0) {
239 maxnumhits1 = maxnumhits0;
242 maxnumhits0 = vtimemap[iclu].size();
244 else if (vtimemap[iclu].
size() > maxnumhits1) {
246 maxnumhits1 = vtimemap[iclu].size();
253 for (
int iclu = 0; iclu < (
int)vtimemap.size(); ++iclu) {
254 if (iclu != iclu1 && iclu != iclu2) iclu3 = iclu;
257 if (iclu1 != -1 && iclu2 != -1) {
259 auto ihit = vhitmap[iclu1].begin();
260 auto itime = vtimemap[iclu1].begin();
261 while (itime != vtimemap[iclu1].
end()) {
264 vtimemap[iclu1].erase(itime++);
265 vhitmap[iclu1].erase(ihit++);
273 ihit = vhitmap[iclu2].begin();
274 itime = vtimemap[iclu2].begin();
275 while (itime != vtimemap[iclu2].
end()) {
278 vtimemap[iclu2].erase(itime++);
279 vhitmap[iclu2].erase(ihit++);
288 if (!vtimemap[iclu1].
size()) {
289 if (iclu3 != -1) {
std::swap(iclu3, iclu1); }
291 if (!vtimemap[iclu2].
size()) {
297 if ((!vtimemap[iclu1].
size()) || (!vtimemap[iclu2].
size()))
return;
300 auto times1 = vtimemap[iclu1].begin();
301 auto timee1 = vtimemap[iclu1].end();
303 auto times2 = vtimemap[iclu2].begin();
304 auto timee2 = vtimemap[iclu2].end();
307 double ts1 = times1->second;
308 double te1 = timee1->second;
309 double ts2 = times2->second;
310 double te2 = timee2->second;
318 double Efield_drift = detProp.Efield(0);
319 double Temperature = detProp.Temperature();
321 double driftvelocity = detProp.DriftVelocity(
322 Efield_drift, Temperature);
323 double timepitch = driftvelocity * timetick;
327 vhitmap[0].begin()->second->WireID().TPC,
328 vhitmap[0].begin()->second->WireID().Cryostat);
331 std::vector<double> vtracklength;
332 for (
size_t iclu = 0; iclu < vtimemap.size(); ++iclu) {
334 double tracklength = 0;
335 if (vtimemap[iclu].
size() == 1) { tracklength = wire_pitch; }
336 else if (!vtimemap[iclu].
empty()) {
338 wend = vtimemap[iclu].cend();
339 double t0 = iw->first, w0 = iw->second;
340 while (++iw != wend) {
341 tracklength +=
std::hypot((iw->first - w0) * wire_pitch, (iw->second - t0) * timepitch);
346 vtracklength.push_back(tracklength);
349 std::map<int, int> maxhitsMatch;
351 auto ihit1 = vhitmap[iclu1].begin();
352 for (
auto itime1 = vtimemap[iclu1].
begin(); itime1 != vtimemap[iclu1].end();
354 std::vector<art::Ptr<recob::Hit>> sp_hits;
355 sp_hits.push_back(ihit1->second);
356 double hitcoord[3] = {0., 0., 0.};
358 if (vtimemap[iclu1].
size() == 1) { length1 = wire_pitch; }
360 for (
auto iw1 = vtimemap[iclu1].
begin(); iw1 != itime1; ++iw1) {
363 length1 +=
std::hypot((iw1->first - iw2->first) * wire_pitch,
364 (iw1->second - iw2->second) * timepitch);
367 double difference = 1e10;
368 auto matchedtime = vtimemap[iclu2].end();
369 auto matchedhit = vhitmap[iclu2].end();
371 auto ihit2 = vhitmap[iclu2].begin();
372 for (
auto itime2 = vtimemap[iclu2].
begin(); itime2 != vtimemap[iclu2].end();
374 if (maxhitsMatch[itime2->first])
continue;
376 if (vtimemap[iclu2].
size() == 1) { length2 = wire_pitch; }
378 for (
auto iw1 = vtimemap[iclu2].
begin(); iw1 != itime2; ++iw1) {
381 length2 +=
std::hypot((iw1->first - iw2->first) * wire_pitch,
382 (iw1->second - iw2->second) * timepitch);
385 if (rev) length2 = vtracklength[iclu2] - length2;
386 length2 = vtracklength[iclu1] / vtracklength[iclu2] * length2;
388 if (timematch &&
std::abs(length2 - length1) < difference) {
389 difference =
std::abs(length2 - length1);
390 matchedtime = itime2;
395 hitcoord[0] = detProp.ConvertTicksToX(
396 matchedtime->second + detProp.GetXTicksOffset((matchedhit->second)->WireID().Plane,
397 (matchedhit->second)->
WireID().TPC,
398 (matchedhit->second)->
WireID().Cryostat),
399 (matchedhit->second)->WireID().Plane,
400 (matchedhit->second)->
WireID().TPC,
401 (matchedhit->second)->
WireID().Cryostat);
415 hitcoord[1] = tmpWIDI.
y;
416 hitcoord[2] = tmpWIDI.
z;
419 if (hitcoord[1] > -1e9 && hitcoord[2] > -1e9) {
420 maxhitsMatch[matchedtime->first] = 1;
421 sp_hits.push_back(matchedhit->second);
424 if (sp_hits.size() > 1) {
425 trajPos.push_back(TVector3(hitcoord));
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
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)
std::vector< TVector3 > trajPos
WireID_t Wire
Index of the wire within its plane.
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void swap(Handle< T > &a, Handle< T > &b)
static int max(int a, int b)
PlaneID_t Plane
Index of the plane within its TPC.
double fsmatch
tolerance for distance matching (in cm)
art::ServiceHandle< geo::Geometry const > geom
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, geo::Point_t &intersection) const
Computes the intersection between two wires.
const detinfo::LArProperties * larprop
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
double y
y position of intersection
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
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
std::vector< std::vector< art::Ptr< recob::Hit > > > trajHit