29 <<
"Passed " << tracks.
size() <<
" tracks for tagging cosmics.";
41 mf::LogInfo(
"pma::PMAlgCosmicTagger") <<
"...tagged " << n <<
" cosmic-like tracks.";
46 mf::LogInfo(
"pma::PMAlgCosmicTagger") <<
" - tag tracks out of 1 drift window;";
49 auto const* geom = lar::providerFrom<geo::Geometry>();
51 for (
auto& t : tracks.
tracks()) {
56 bool node_out_of_drift_min =
false;
57 bool node_out_of_drift_max =
false;
58 for (
size_t nidx = 0; nidx < trk.
Nodes().size(); ++nidx) {
59 auto const& node = *(trk.
Nodes()[nidx]);
60 auto const& tpcGeo = geom->TPC(
geo::TPCID(node.Cryo(), node.TPC()));
62 int driftDir =
abs(tpcGeo.DetectDriftDirection());
63 p = node.Point3D()[driftDir - 1];
79 <<
"Drift direction unknown: " << driftDir << std::endl;
85 if (node_out_of_drift_min && node_out_of_drift_max) {
90 else if (node_out_of_drift_min || node_out_of_drift_max) {
97 mf::LogInfo(
"pma::PMAlgCosmicTagger") <<
" - Tagged " << n <<
" tracks out of 1 drift window.";
110 for (
auto& t : tracks.
tracks()) {
113 if (t.Track()->GetT0() != 0.0) {
114 mf::LogInfo(
"pma::PMAlgCosmicTagger") <<
" - track with T0 = " << t.Track()->GetT0();
124 mf::LogInfo(
"pma::PMAlgCosmicTagger") <<
" - Tagged " << n <<
" non-beam T0 tracks.";
133 auto const* geom = lar::providerFrom<geo::Geometry>();
139 for (
auto& t : tracks.
tracks()) {
142 auto const& node0 = *(t.Track()->Nodes()[0]);
143 auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);
147 (node0.Point3D()[hIdx] > node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
149 (node0.Point3D()[hIdx] <= node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
158 if (top && frontBack) {
166 <<
" - Tagged " << n <<
" tracks crossing from top to front/back." << std::endl;
179 auto const* geom = lar::providerFrom<geo::Geometry>();
184 for (
auto& t : tracks.
tracks()) {
187 auto const& node0 = *(t.Track()->Nodes()[0]);
188 auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);
192 (node0.Point3D()[hIdx] > node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
194 (node0.Point3D()[hIdx] <= node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
201 bool foundTrack =
false;
202 for (
auto const&
tt : tracks.
tracks()) {
204 if ((&
tt) == (&t))
continue;
207 TVector3 trkVtx = (
tt.Track()->Nodes()[0])->
Point3D();
208 TVector3 trkEnd = (
tt.Track()->Nodes()[
tt.Track()->Nodes().size() - 1])->
Point3D();
230 std::vector<float> nSigmaPerView;
233 for (
auto const view : geom->Views()) {
236 std::map<size_t, std::vector<double>> dedx;
237 t.Track()->GetRawdEdxSequence(dedx, view);
239 std::vector<double> trk_dedx;
241 for (
int h = t.Track()->NextHit(-1, view); h != -1; h = t.Track()->NextHit(h, view)) {
243 if (h > t.Track()->PrevHit(t.Track()->size(), view))
break;
245 if (dedx[h][5] / dedx[h][6] <= 0 || dedx[h][5] / dedx[h][6] > 1e6)
continue;
246 trk_dedx.push_back(dedx[h][5] / dedx[h][6]);
249 if (trk_dedx.size() == 0) {
251 <<
"View " << view <<
" has no hits." << std::endl;
256 double mean = sum /
static_cast<double>(trk_dedx.size());
259 accum += (d -
mean) * (d - mean);
261 double stdev = sqrt(accum / static_cast<double>(trk_dedx.size() - 1));
264 <<
" View " << view <<
" has average dedx " << mean <<
" +/- " << stdev
265 <<
" and final dedx " << trk_dedx[trk_dedx.size() - 1] << std::endl;
267 nSigmaPerView.push_back(fabs((trk_dedx[trk_dedx.size() - 1] -
mean) / stdev));
270 bool notStopper =
true;
271 short unsigned int n2Sigma = 0;
272 short unsigned int n3Sigma = 0;
273 for (
auto const nSigma : nSigmaPerView) {
274 if (nSigma >= 2.0) ++n2Sigma;
275 if (nSigma >= 3.0) ++n3Sigma;
278 if (n3Sigma > 0) notStopper =
false;
279 if (n2Sigma == nSigmaPerView.size()) notStopper =
false;
291 <<
" - Tagged " << n <<
" tracks stopping in the detector after starting at the top." 301 auto const* geom = lar::providerFrom<geo::Geometry>();
302 auto const dir = geom->TPC().HeightDir();
307 <<
" - Tagged " << n <<
" tracks crossing the full detector height";
315 auto const* geom = lar::providerFrom<geo::Geometry>();
316 auto const dir = geom->TPC().WidthDir();
321 <<
" - Tagged " << n <<
" tracks crossing the full detector width";
329 auto const* geom = lar::providerFrom<geo::Geometry>();
330 auto const dir = geom->TPC().LengthDir();
335 <<
" - Tagged " << n <<
" tracks crossing the full detector length";
343 if (direction == -1) {
345 <<
" - Could not recognise direction, not attempting to perform fullCrossingTagger.";
362 for (
auto& t : tracks.
tracks()) {
365 auto const& node0 = *(t.Track()->Nodes()[0]);
366 auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);
369 double trkDim = fabs(node0.Point3D()[direction] - node1.Point3D()[direction]);
374 t.Track()->SetTagFlag(dirTag);
375 mf::LogInfo(
"pma::PMAlgCosmicTagger") <<
" -- track tagged in direction " << direction
376 <<
" with " << trkDim <<
" (c.f. " << detDim <<
")";
385 short int dirIndx)
const 388 return (fabs(pos[dirIndx] -
fDimensionsMax[dirIndx]) < tolerance);
393 short int dirIndx)
const 396 bool front = (fabs(pos[dirIndx] -
fDimensionsMin[dirIndx]) < tolerance);
397 bool back = (fabs(pos[dirIndx] -
fDimensionsMax[dirIndx]) < tolerance);
399 return front || back;
413 auto const* geom = lar::providerFrom<geo::Geometry>();
420 if (TPC.DriftDistance() < 25.0) {
continue; }
423 auto const center = TPC.GetCenter();
424 double tpcDim[3] = {TPC.HalfWidth(), TPC.HalfHeight(), 0.5 * TPC.Length()};
426 if (center.X() - tpcDim[0] < minX) minX = center.X() - tpcDim[0];
427 if (center.X() + tpcDim[0] > maxX) maxX = center.X() + tpcDim[0];
428 if (center.Y() - tpcDim[1] <
minY) minY = center.Y() - tpcDim[1];
429 if (center.Y() + tpcDim[1] >
maxY) maxY = center.Y() + tpcDim[1];
430 if (center.Z() - tpcDim[2] < minZ) minZ = center.Z() - tpcDim[2];
431 if (center.Z() + tpcDim[2] > maxZ) maxZ = center.Z() + tpcDim[2];
447 if (dir.X() > 0.99)
return 0;
448 if (dir.Y() > 0.99)
return 1;
size_t fullWidthCrossing(pma::TrkCandidateColl &tracks) const
Utilities related to art service access.
bool fTagOutOfDriftTracks
void SetTagFlag(ETag value)
Implementation of the Projection Matching Algorithm.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Geometry information for a single TPC.
::fhicl::TupleAs< Point(::geo::Length_t,::geo::Length_t,::geo::Length_t)> Point3D
Atom object for reading a 3D point or vector (centimeters).
double fTopFrontBackMargin
bool fTagFullLengthTracks
constexpr auto abs(T v)
Returns the absolute value of the argument.
Implementation of the Projection Matching Algorithm.
short int ConvertDirToInt(const geo::Vector_t &dir) const
std::vector< double > fDimensionsMin
size_t tagApparentStopper(pma::TrkCandidateColl &tracks) const
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Access the description of detector geometry.
size_t outOfDriftWindow(pma::TrkCandidateColl &tracks) const
bool isFrontBackVertex(const TVector3 &pos, double tolerance, short int dirIndx) const
bool fTagFullHeightTracks
size_t nonBeamT0Tag(detinfo::DetectorClocksData const &clockData, pma::TrkCandidateColl &tracks) const
double fFullCrossingMargin
std::vector< double > fDimensionsMax
The data type to uniquely identify a TPC.
void tag(detinfo::DetectorClocksData const &clockData, pma::TrkCandidateColl &tracks)
Track finding helper for the Projection Matching Algorithm.
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Contains all timing reference information for the detector.
size_t tagTopFrontBack(pma::TrkCandidateColl &tracks) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
int trigger_offset(DetectorClocksData const &data)
size_t fullHeightCrossing(pma::TrkCandidateColl &tracks) const
size_t fullCrossingTagger(pma::TrkCandidateColl &tracks, int direction) const
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
std::vector< pma::Node3D * > const & Nodes() const noexcept
bool isTopVertex(const TVector3 &pos, double tolerance, short int dirIndx) const
size_t fullLengthCrossing(pma::TrkCandidateColl &tracks) const
double fApparentStopperMargin
art framework interface to geometry description
std::vector< TrkCandidate > const & tracks() const
cet::coded_exception< error, detail::translate > exception
Encapsulate the construction of a single detector plane.