26 std::vector<TVector3>& TrajPos,
27 std::vector<TVector3>& TrajDir)
59 unsigned short minLen = 9999;
60 unsigned short maxLen = 0;
61 unsigned short nPlnsWithHits = 0;
62 unsigned short ipl, aPlane = 3;
63 for (ipl = 0; ipl < 3; ++ipl) {
64 if (trkX[ipl].
size() == 0)
continue;
66 if (trkX[ipl].
size() < minLen) { minLen = trkX[ipl].size(); }
67 if (trkX[ipl].
size() > maxLen) {
68 maxLen = trkX[ipl].size();
73 mf::LogVerbatim(
"TTA") <<
"trkX sizes " << trkX[0].size() <<
" " << trkX[1].size() <<
" " 74 << trkX[2].size() <<
" " 75 <<
" nPlnsWithHits " << nPlnsWithHits;
76 if (nPlnsWithHits < 2)
return;
77 if (aPlane > 2)
return;
87 iht = trkX[aPlane].size() - 1;
88 if (trkX[aPlane][0] > trkX[aPlane][iht]) {
90 for (ipl = 0; ipl < 3; ++ipl) {
91 if (trkX[ipl].
size() == 0)
continue;
92 std::reverse(trkWID[ipl].
begin(), trkWID[ipl].
end());
93 std::reverse(trkX[ipl].
begin(), trkX[ipl].
end());
94 std::reverse(trkXErr[ipl].
begin(), trkXErr[ipl].
end());
101 myprt <<
"TrkXW end0";
102 for (ipl = 0; ipl < 3; ++ipl) {
103 if (trkX[ipl].
size() == 0)
continue;
104 myprt <<
" " << trkX[ipl][0];
107 myprt <<
"TrkXW end1";
108 for (ipl = 0; ipl < 3; ++ipl) {
109 if (trkX[ipl].
size() == 0)
continue;
110 iht = trkX[ipl].size() - 1;
111 myprt <<
" " << trkX[ipl][iht];
122 for (
unsigned short ipl = 0; ipl < 3; ++ipl) {
123 if (trkX[ipl].
size() == 0)
continue;
124 if (trkX[ipl][0] <
minX) {
128 iht = trkX[ipl].size() - 1;
129 if (trkX[ipl][iht] >
maxX) {
130 maxX = trkX[ipl][iht];
141 double aveHitXErr = 0;
142 unsigned short nHit = 0;
143 for (ipl = 0; ipl < 3; ++ipl) {
144 for (iht = 0; iht < trkXErr[ipl].size(); ++iht) {
145 aveHitXErr += trkXErr[ipl][iht];
149 aveHitXErr /= (double)nHit;
151 unsigned short npt = (
maxX -
minX) / (1 * aveHitXErr);
152 if (npt < 2) npt = 2;
153 if (npt > maxLen) npt = maxLen;
155 if (
prt)
mf::LogVerbatim(
"TTA") <<
" aveHitXErr " << aveHitXErr <<
" number of traj points ";
157 double maxBinX = (
maxX -
minX) / (
double)(npt - 1);
158 double binX = maxBinX;
159 double xOrigin =
minX;
163 std::vector<geo::WireID> hitWID;
164 std::vector<double> hitX;
165 std::vector<double> hitXErr;
169 std::vector<TVector3> STPos;
170 std::vector<TVector3> STDir;
173 if (STPos.size() != 2 || STDir.size() != STPos.size()) {
177 for (ipl = 0; ipl < 3; ++ipl) {
178 if (trkX[ipl].
size() == 0)
continue;
179 std::reverse(trkWID[ipl].
begin(), trkWID[ipl].
end());
180 std::reverse(trkX[ipl].
begin(), trkX[ipl].
end());
181 std::reverse(trkXErr[ipl].
begin(), trkXErr[ipl].
end());
189 for (
unsigned short ii = 0; ii < STPos.size(); ++ii)
190 mf::LogVerbatim(
"TTA") << ii <<
" " << std::fixed << std::setprecision(1) << STPos[ii](0)
191 <<
" " << STPos[ii](1) <<
" " << STPos[ii](2);
194 if (maxLen < 4 || npt < 2) {
199 std::reverse(TrajPos.begin(), TrajPos.end());
200 std::reverse(TrajDir.begin(), TrajDir.end());
201 for (ipl = 0; ipl < 3; ++ipl) {
202 if (trkX[ipl].
size() == 0)
continue;
203 std::reverse(trkWID[ipl].
begin(), trkWID[ipl].
end());
204 std::reverse(trkX[ipl].
begin(), trkX[ipl].
end());
205 std::reverse(trkXErr[ipl].
begin(), trkXErr[ipl].
end());
212 std::array<unsigned short, 3> hStart;
213 for (ipl = 0; ipl < 3; ++ipl)
216 bool gotLastPoint =
false;
217 for (
unsigned short ipt = 0; ipt < npt + 1; ++ipt) {
221 for (ipl = 0; ipl < 3; ++ipl) {
222 for (iht = hStart[ipl]; iht < trkX[ipl].size(); ++iht) {
223 if (trkX[ipl][iht] < xOrigin - binX)
continue;
224 if (trkX[ipl][iht] > xOrigin + binX)
break;
225 hitWID.push_back(trkWID[ipl][iht]);
226 hitX.push_back(trkX[ipl][iht]);
227 hitXErr.push_back(trkXErr[ipl][iht]);
232 mf::LogVerbatim(
"TTA") <<
"ipt " << ipt <<
" xOrigin " << xOrigin <<
" binX " << binX
233 <<
" hitX size " << hitX.size();
234 if (hitX.size() > 3) {
237 mf::LogVerbatim(
"TTA") <<
" xyz " << xyz(0) <<
" " << xyz(1) <<
" " << xyz(2) <<
" dir " 238 <<
dir(0) <<
" " <<
dir(1) <<
" " <<
dir(2) <<
" ChiDOF " << ChiDOF
239 <<
" hitX size " << hitX.size();
241 else if (ipt == 0 && STPos.size() == 2) {
247 if (xOrigin >=
maxX)
break;
249 if (
maxX - xOrigin < binX) { xOrigin =
maxX; }
253 if (ChiDOF < 0 || ChiDOF > 100)
continue;
254 TrajPos.push_back(xyz);
255 TrajDir.push_back(dir);
256 if (ipt == npt) gotLastPoint =
true;
259 mf::LogVerbatim(
"TTA") <<
"gotLastPoint " << gotLastPoint <<
" TTA Traj \n";
260 for (
unsigned short ii = 0; ii < TrajPos.size(); ++ii)
261 mf::LogVerbatim(
"TTA") << ii <<
" " << std::fixed << std::setprecision(1) << TrajPos[ii](0)
262 <<
" " << TrajPos[ii](1) <<
" " << TrajPos[ii](2);
265 if (!gotLastPoint && STPos.size() == 2) {
267 TrajPos.push_back(STPos[1]);
268 TrajDir.push_back(STDir[1]);
271 if (TrajPos.size() < 2) {
278 std::reverse(TrajPos.begin(), TrajPos.end());
279 std::reverse(TrajDir.begin(), TrajDir.end());
280 for (ipl = 0; ipl < 3; ++ipl) {
281 if (trkX[ipl].
size() == 0)
continue;
282 std::reverse(trkWID[ipl].
begin(), trkWID[ipl].
end());
283 std::reverse(trkX[ipl].
begin(), trkX[ipl].
end());
284 std::reverse(trkXErr[ipl].
begin(), trkXErr[ipl].
end());
290 for (
unsigned short ii = 0; ii < TrajPos.size(); ++ii)
291 mf::LogVerbatim(
"TTA") << ii <<
" " << std::fixed << std::setprecision(1) << TrajPos[ii](0)
292 <<
" " << TrajPos[ii](1) <<
" " << TrajPos[ii](2);
301 std::vector<TVector3>& TrajPos,
302 std::vector<TVector3>& TrajDir)
306 std::vector<unsigned short> usePlns;
307 double y,
z, endX, endY, endZ;
308 unsigned short ipl, iht, nend, jpl, iPlane, jPlane;
309 unsigned int iWire, jWire;
313 for (
unsigned short nit = 0; nit < 5; ++nit) {
315 for (ipl = 0; ipl < 3; ++ipl) {
316 if (trkX[ipl].
size() == 0)
continue;
317 if (trkX[ipl][0] < xCut) usePlns.push_back(ipl);
319 if (usePlns.size() >= 2)
break;
324 mf::LogVerbatim(
"TTA") <<
"ShortTrack minX end " << xCut <<
" usePlns size " 326 if (usePlns.size() < 2)
return;
332 endX = trkX[ipl][iht];
333 iPlane = trkWID[ipl][iht].Plane;
334 iWire = trkWID[ipl][iht].Wire;
335 unsigned int cstat = trkWID[ipl][iht].Cryostat;
336 unsigned int tpc = trkWID[ipl][iht].TPC;
337 for (
unsigned short jj = 1; jj < usePlns.size(); ++jj) {
339 endX += trkX[jpl][iht];
340 jPlane = trkWID[jpl][iht].Plane;
341 jWire = trkWID[jpl][iht].Wire;
350 xyz(0) = endX / (float)(nend + 1);
351 xyz(1) = endY / (float)nend;
352 xyz(2) = endZ / (float)nend;
353 if (
prt)
mf::LogVerbatim(
"TTA") <<
" xyz " << xyz(0) <<
" " << xyz(1) <<
" " << xyz(2);
354 TrajPos.push_back(xyz);
360 for (
unsigned short nit = 0; nit < 5; ++nit) {
362 for (ipl = 0; ipl < 3; ++ipl) {
363 if (trkX[ipl].
size() == 0)
continue;
364 iht = trkX[ipl].size() - 1;
365 if (trkX[ipl][iht] > xCut) usePlns.push_back(ipl);
367 if (usePlns.size() >= 2)
break;
372 mf::LogVerbatim(
"TTA") <<
"ShortTrack maxX end " << xCut <<
" usePlns size " 374 if (usePlns.size() < 2) {
383 iht = trkX[ipl].size() - 1;
384 endX = trkX[ipl][iht];
385 iPlane = trkWID[ipl][iht].Plane;
386 iWire = trkWID[ipl][iht].Wire;
387 for (
unsigned short jj = 1; jj < usePlns.size(); ++jj) {
389 iht = trkX[jpl].size() - 1;
390 endX += trkX[jpl][iht];
391 jPlane = trkWID[jpl][iht].Plane;
392 jWire = trkWID[jpl][iht].Wire;
400 xyz(0) = endX / (float)(nend + 1);
401 xyz(1) = endY / (float)nend;
402 xyz(2) = endZ / (float)nend;
403 if (
prt)
mf::LogVerbatim(
"TTA") <<
" xyz " << xyz(0) <<
" " << xyz(1) <<
" " << xyz(2);
404 TrajPos.push_back(xyz);
405 TVector3
dir = TrajPos[1] - TrajPos[0];
407 TrajDir.push_back(dir);
408 TrajDir.push_back(dir);
411 mf::LogVerbatim(
"TTA") <<
">>>> Short track (" << TrajPos[0](0) <<
", " << TrajPos[0](1)
412 <<
", " << TrajPos[0](2) <<
") to (" << TrajPos[1](0) <<
", " 413 << TrajPos[1](1) <<
", " << TrajPos[1](2) <<
")";
void TrkLineFit(std::vector< geo::WireID > &hitWID, std::vector< double > &hitX, std::vector< double > &hitXErr, double XOrigin, TVector3 &Pos, TVector3 &Dir, float &ChiDOF) const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
TrackLineFitAlg fTrackLineFitAlg
unsigned short fMaxTrajPoints
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 array(Array const &a)
Returns a manipulator which will print the specified array.
void TrackTrajectory(std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
Definition of data types for geometry description.
void ShortTrackTrajectory(std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
art::ServiceHandle< geo::Geometry const > geom
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.