32 bool greaterThan1(PlnLen p1, PlnLen p2)
34 return p1.length > p2.length;
42 fSPTAlg = pset.
get<
int>(
"SPTAlg", 0);
43 fTrajOnly = pset.
get<
bool>(
"TrajOnly",
false);
44 ftmatch = pset.
get<
double>(
"TMatch");
45 fsmatch = pset.
get<
double>(
"SMatch");
60 if (fSPTAlg == 0) { TrackTrajectory(detProp, fHits); }
61 else if (fSPTAlg == 1) {
62 Track3D(clockData, detProp, fHits);
66 <<
"Unknown SPTAlg " << fSPTAlg <<
", needs to be 0 or 1";
69 if (!fTrajOnly && trajPos.size())
70 MakeSPT(clockData, detProp, fHits);
82 std::array<std::vector<geo::WireID>, 3> trkWID;
83 std::array<std::vector<double>, 3> trkX;
84 std::array<std::vector<double>, 3> trkXErr;
86 std::array<std::vector<art::Ptr<recob::Hit>>, 3> trajHits;
88 for (
size_t i = 0; i < fHits.size(); ++i) {
89 trajHits[fHits[i]->WireID().Plane].push_back(fHits[i]);
92 std::vector<PlnLen> spl;
94 for (
unsigned int ipl = 0; ipl < 3; ++ipl) {
96 plnlen.length = trajHits[ipl].size();
97 spl.push_back(plnlen);
99 std::sort(spl.begin(), spl.end(), greaterThan1);
102 for (
size_t ipl = 0; ipl < 3; ++ipl) {
103 if (!trajHits[ipl].
size())
continue;
104 double xbeg0 = detProp.
ConvertTicksToX(trajHits[spl[0].index][0]->PeakTime(),
106 trajHits[spl[0].index][0]->
WireID().TPC,
107 trajHits[spl[0].index][0]->
WireID().Cryostat);
108 double xend0 = detProp.
ConvertTicksToX(trajHits[spl[0].index].back()->PeakTime(),
110 trajHits[spl[0].index].back()->
WireID().TPC,
111 trajHits[spl[0].index].back()->
WireID().Cryostat);
114 trajHits[ipl][0]->
WireID().TPC,
115 trajHits[ipl][0]->
WireID().Cryostat);
116 double xend1 = detProp.
ConvertTicksToX(trajHits[ipl].back()->PeakTime(),
118 trajHits[ipl].back()->
WireID().TPC,
119 trajHits[ipl].back()->
WireID().Cryostat);
126 std::reverse(trajHits[ipl].
begin(), trajHits[ipl].
end());
130 for (
size_t ipl = 0; ipl < 3; ++ipl) {
132 trkWID[ipl].resize(trajHits[ipl].
size());
133 trkX[ipl].resize(trajHits[ipl].
size());
134 trkXErr[ipl].resize(trajHits[ipl].
size());
135 for (
size_t iht = 0; iht < trajHits[ipl].size(); ++iht) {
136 trkWID[ipl][iht] = trajHits[ipl][iht]->WireID();
137 trkX[ipl][iht] = detProp.
ConvertTicksToX(trajHits[ipl][iht]->PeakTime(),
139 trajHits[ipl][iht]->
WireID().TPC,
140 trajHits[ipl][iht]->
WireID().Cryostat);
141 trkXErr[ipl][iht] = 0.2 * trajHits[ipl][iht]->RMS() * trajHits[ipl][iht]->Multiplicity();
144 fTrackTrajectoryAlg.TrackTrajectory(trkWID, trkX, trkXErr, trajPos, trajDir);
145 if (!trajPos.size() || !trajDir.size()) {
146 mf::LogWarning(
"CosmicTrackerAlg") <<
"Failed to reconstruct trajectory points.";
150 for (
auto ipos = trajPos.begin(), idir = trajDir.begin(); ipos != trajPos.end() - 1;) {
151 auto ipos1 = ipos + 1;
152 if (*ipos1 == *ipos) {
153 ipos = trajPos.erase(ipos);
154 idir = trajDir.erase(idir);
164 vw.resize(geom->Ncryostats());
165 vt.resize(geom->Ncryostats());
166 vtraj.resize(geom->Ncryostats());
168 auto const cstat = cryostat.ID().Cryostat;
169 vw[cstat].resize(cryostat.NTPC());
170 vt[cstat].resize(cryostat.NTPC());
171 vtraj[cstat].resize(cryostat.NTPC());
172 for (
size_t tpc = 0; tpc < cryostat.NTPC(); ++tpc) {
174 vw[cstat][tpc].resize(tpcgeom.
Nplanes());
175 vt[cstat][tpc].resize(tpcgeom.
Nplanes());
176 vtraj[cstat][tpc].resize(tpcgeom.
Nplanes());
177 for (
unsigned int plane = 0; plane < tpcgeom.
Nplanes(); ++plane) {
178 for (
size_t i = 0; i < trajPos.size(); ++i) {
179 double wirecord = geom->WireCoordinate(
geo::Point_t{0, trajPos[i].Y(), trajPos[i].Z()},
182 vw[cstat][tpc][plane].push_back(wirecord);
183 vt[cstat][tpc][plane].push_back(tick);
184 vtraj[cstat][tpc][plane].push_back(i);
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();
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();
247 std::swap(iclu1, iclu2);
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()) {
290 std::swap(iclu3, iclu2);
291 std::swap(iclu1, iclu2);
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);
319 Efield_drift, Temperature);
320 double timepitch = driftvelocity * timetick;
323 geom->WirePitch(vhitmap[0].
begin()->
second->WireID().asPlaneID());
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;
382 bool timematch =
std::abs(itime1->second - itime2->second) < ftmatch;
383 if (timematch &&
std::abs(length2 - length1) < difference) {
384 difference =
std::abs(length2 - length1);
385 matchedtime = itime2;
389 if (difference < fsmatch) {
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);
407 bool sameTpcOrNot = geom->WireIDsIntersect(c1, c2, tmpWIDI);
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));
421 trajHit.push_back(sp_hits);
432 larprop = lar::providerFrom<detinfo::LArPropertiesService>();
435 double Efield_drift = detProp.
Efield(0);
437 double driftvelocity =
439 double time_pitch = driftvelocity * timetick;
441 for (
size_t i = 0; i < fHits.size(); ++i) {
448 double mindis1 = 1e9;
449 double mindis2 = 1e9;
450 double mindis2p = 1e9;
452 unsigned int wire = wireid.
Wire;
453 unsigned int plane = wireid.
Plane;
454 unsigned int tpc = wireid.
TPC;
455 unsigned int cstat = wireid.
Cryostat;
456 double wire_pitch = geom->WirePitch(wireid.
asPlaneID());
459 for (
size_t j = 0; j < vw[cstat][tpc][plane].size(); ++j) {
460 double dis = std::hypot(wire_pitch * (wire - vw[cstat][tpc][plane][j]),
461 time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]));
464 ip1 = vtraj[cstat][tpc][plane][j];
469 for (
size_t j = 0; j < vw[cstat][tpc][plane].size(); ++j) {
470 if (
int(j) == ih1)
continue;
471 double dis = std::hypot(wire_pitch * (wire - vw[cstat][tpc][plane][j]),
472 time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]));
475 ip2 = vtraj[cstat][tpc][plane][j];
478 TVector3 v1(wire_pitch * (wire - vw[cstat][tpc][plane][ih1]),
479 time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]),
481 TVector3 v2(wire_pitch * (wire - vw[cstat][tpc][plane][j]),
482 time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][j]),
484 if (v1.Dot(v2) < 0) {
485 if (dis < mindis2p) {
487 ip2p = vtraj[cstat][tpc][plane][j];
496 if (ip1 == -1 || ip2 == -1) {
498 throw cet::exception(
"CosmicTrackerAlg") <<
"Cannot find two nearest trajectory points.";
500 TVector3 v1(wire_pitch * (wire - vw[cstat][tpc][plane][ih1]),
501 time_pitch * (fHits[i]->PeakTime() - vt[cstat][tpc][plane][ih1]),
503 TVector3 v2(wire_pitch * (vw[cstat][tpc][plane][ih2] - vw[cstat][tpc][plane][ih1]),
504 time_pitch * (vt[cstat][tpc][plane][ih2] - vt[cstat][tpc][plane][ih1]),
507 throw cet::exception(
"CosmicTrackerAlg") <<
"Two identical trajectory points.";
509 double ratio = v1.Dot(v2) / v2.Mag2();
510 TVector3 pos = trajPos[ip1] + ratio * (trajPos[ip2] - trajPos[ip1]);
511 trkPos.push_back(pos);
512 if (trajDir[ip1].Mag()) { trkDir.push_back(trajDir[ip1]); }
513 else if (trajDir[ip2].Mag()) {
514 trkDir.push_back(trajDir[ip2]);
517 trkDir.push_back((trajPos[ip2] - trajPos[ip1]).Unit());
geo::TPCID const & ID() const
Returns the identifier of this TPC.
code to link reconstructed objects back to the MC truth information
double z
z position of intersection
Utilities related to art service access.
Encapsulate the construction of a single cyostat.
double GetXTicksOffset(int p, int t, int c) const
unsigned int Nplanes() const
Number of planes in this tpc.
Declaration of signal hit object.
The data type to uniquely identify a Plane.
Geometry information for a single TPC.
double Temperature() const
In kelvin.
constexpr auto abs(T v)
Returns the absolute value of the argument.
CryostatID_t Cryostat
Index of cryostat.
WireID_t Wire
Index of the wire within its plane.
Geometry information for a single cryostat.
double Efield(unsigned int planegap=0) const
kV/cm
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
T get(std::string const &key) const
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
double ConvertXToTicks(double X, int p, int t, int c) const
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
PlaneID_t Plane
Index of the plane within its TPC.
Definition of data types for geometry description.
Class defining a plane for tracking. It provides various functionalities to convert track parameters ...
double ConvertTicksToX(double ticks, int p, int t, int c) const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Contains all timing reference information for the detector.
double y
y position of intersection
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
constexpr WireID()=default
Default constructor: an invalid TPC ID.
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
TPCID_t TPC
Index of the TPC within its cryostat.
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.
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Encapsulate the construction of a single detector plane.