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();
250 std::swap(iclu1, iclu2);
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()) {
262 if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) -
ftmatch ||
263 itime->second > std::min(tmax[iclu1], tmax[iclu2]) +
ftmatch) {
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()) {
276 if (itime->second < std::max(tmin[iclu1], tmin[iclu2]) -
ftmatch ||
277 itime->second > std::min(tmax[iclu1], tmax[iclu2]) +
ftmatch) {
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()) {
293 std::swap(iclu3, iclu2);
294 std::swap(iclu1, iclu2);
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()) {
337 std::map<int, double>::const_iterator iw = vtimemap[iclu].cbegin(),
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 = 1
e10;
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);
412 bool sameTpcOrNot =
geom->WireIDsIntersect(c1, c2, tmpWIDI);
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));
double z
z position of intersection
double ftmatch
tolerance for time matching (in ticks)
std::vector< TVector3 > trajPos
std::size_t size(FixedBins< T, C > const &) noexcept
WireID_t Wire
Index of the wire within its plane.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
auto end(FixedBins< T, C > const &) noexcept
PlaneID_t Plane
Index of the plane within its TPC.
double fsmatch
tolerance for distance matching (in cm)
auto begin(FixedBins< T, C > const &) noexcept
art::ServiceHandle< geo::Geometry const > geom
const detinfo::LArProperties * larprop
double y
y position of intersection
then echo ***************************************echo Variable FHICL_FILE_PATH not found echo You porbably haven t set up larsoft echo Try setup uboonecode vXX_XX_XX q e10
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
bool empty(FixedBins< T, C > const &) noexcept
std::vector< std::vector< art::Ptr< recob::Hit > > > trajHit