pma::PMAlgCosmicTagger Class Reference

#include <PMAlgCosmicTagger.h>

Classes
struct	Config

Public Member Functions
	PMAlgCosmicTagger (const pma::PMAlgCosmicTagger::Config &config)
bool	tagAny () const
void	tag (detinfo::DetectorClocksData const &clockData, pma::TrkCandidateColl &tracks)

Private Member Functions
size_t	outOfDriftWindow (pma::TrkCandidateColl &tracks) const
size_t	fullHeightCrossing (pma::TrkCandidateColl &tracks) const
size_t	fullWidthCrossing (pma::TrkCandidateColl &tracks) const
size_t	fullLengthCrossing (pma::TrkCandidateColl &tracks) const
size_t	fullCrossingTagger (pma::TrkCandidateColl &tracks, int direction) const
size_t	nonBeamT0Tag (detinfo::DetectorClocksData const &clockData, pma::TrkCandidateColl &tracks) const
size_t	tagTopFrontBack (pma::TrkCandidateColl &tracks) const
size_t	tagApparentStopper (pma::TrkCandidateColl &tracks) const
bool	isTopVertex (const TVector3 &pos, double tolerance, short int dirIndx) const
bool	isFrontBackVertex (const TVector3 &pos, double tolerance, short int dirIndx) const
void	GetDimensions ()
short int	ConvertDirToInt (const TVector3 &dir) const

Private Attributes
bool	fTagOutOfDriftTracks
double	fOutOfDriftMargin
bool	fTagFullHeightTracks
bool	fTagFullWidthTracks
bool	fTagFullLengthTracks
double	fFullCrossingMargin
bool	fTagNonBeamT0Tracks
double	fNonBeamT0Margin
bool	fTagTopFrontBack
double	fTopFrontBackMargin
bool	fTagApparentStopper
double	fApparentStopperMargin
bool	fVetoActualStopper
double	fStopperBuffer
std::vector< double >	fDimensionsMin
std::vector< double >	fDimensionsMax

Detailed Description

Definition at line 27 of file PMAlgCosmicTagger.h.

Constructor & Destructor Documentation

pma::PMAlgCosmicTagger::PMAlgCosmicTagger ( const pma::PMAlgCosmicTagger::Config &  config )

inline

Definition at line 80 of file PMAlgCosmicTagger.h.

    : fTagOutOfDriftTracks(config.TagOutOfDriftTracks())
    , fOutOfDriftMargin(config.OutOfDriftMargin())
    ,

    fTagFullHeightTracks(config.TagFullHeightTracks())
    , fTagFullWidthTracks(config.TagFullWidthTracks())
    , fTagFullLengthTracks(config.TagFullLengthTracks())
    , fFullCrossingMargin(config.FullCrossingMargin())
    ,

    fTagNonBeamT0Tracks(config.TagNonBeamT0Tracks())
    , fNonBeamT0Margin(config.NonBeamT0Margin())
    , fTagTopFrontBack(config.TagTopFrontBack())
    , fTopFrontBackMargin(config.TopFrontBackMargin())
    , fTagApparentStopper(config.TagApparentStopper())
    , fApparentStopperMargin(config.ApparentStopperMargin())
    , fVetoActualStopper(config.VetoActualStopper())
    , fStopperBuffer(config.StopperBuffer())
  {}

Member Function Documentation

short int pma::PMAlgCosmicTagger::ConvertDirToInt ( const TVector3 &  dir ) const

private

Definition at line 458 of file PMAlgCosmicTagger.cxx.

{

  if (dir.X() > 0.99) return 0;
  if (dir.Y() > 0.99) return 1;
  if (dir.Z() > 0.99)
    return 2;

  else
    return -1;
}

size_t pma::PMAlgCosmicTagger::fullCrossingTagger ( pma::TrkCandidateColl &  tracks, int  direction  ) const

private

Definition at line 348 of file PMAlgCosmicTagger.cxx.

{

  if (direction == -1) {
    mf::LogWarning("pma::PMAlgCosmicTagger")
      << " - Could not recognise direction, not attempting to perform fullCrossingTagger.";
    return 0;
  }

  size_t n = 0;

  double detDim = fDimensionsMax[direction] - fDimensionsMin[direction];

  pma::Track3D::ETag dirTag = pma::Track3D::kNotTagged;
  switch (direction) {
  case 0: dirTag = pma::Track3D::kGeometry_XX; break;
  case 1: dirTag = pma::Track3D::kGeometry_YY; break;
  case 2: dirTag = pma::Track3D::kGeometry_ZZ; break;
  default: dirTag = pma::Track3D::kNotTagged; break;
  }

  // Loop over the tracks
  for (auto& t : tracks.tracks()) {

    // Get the first and last positions from the track.
    auto const& node0 = *(t.Track()->Nodes()[0]);
    auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);

    // Get the length of the track in the requested direction
    double trkDim = fabs(node0.Point3D()[direction] - node1.Point3D()[direction]);

    if ((detDim - trkDim) < fFullCrossingMargin) {
      ++n;
      t.Track()->SetTagFlag(pma::Track3D::kCosmic);
      t.Track()->SetTagFlag(dirTag);
      mf::LogInfo("pma::PMAlgCosmicTagger") << " -- track tagged in direction " << direction
                                            << " with " << trkDim << " (c.f. " << detDim << ")";
    }
  }

  return n;
}

size_t pma::PMAlgCosmicTagger::fullHeightCrossing ( pma::TrkCandidateColl &  tracks ) const

private

Definition at line 303 of file PMAlgCosmicTagger.cxx.

{

  // Just use the first tpc to get the height dimension (instead of assuming y).
  auto const* geom = lar::providerFrom<geo::Geometry>();
  TVector3 dir = geom->TPC(0, 0).HeightDir();

  size_t n = fullCrossingTagger(tracks, ConvertDirToInt(dir));

  mf::LogInfo("pma::PMAlgCosmicTagger")
    << " - Tagged " << n << " tracks crossing the full detector height";
  return n;
}

size_t pma::PMAlgCosmicTagger::fullLengthCrossing ( pma::TrkCandidateColl &  tracks ) const

private

Definition at line 333 of file PMAlgCosmicTagger.cxx.

{

  // Just use the first tpc to get the length dimension (instead of assuming z).
  auto const* geom = lar::providerFrom<geo::Geometry>();
  TVector3 dir = geom->TPC(0, 0).LengthDir();

  size_t n = fullCrossingTagger(tracks, ConvertDirToInt(dir));

  mf::LogInfo("pma::PMAlgCosmicTagger")
    << " - Tagged " << n << " tracks crossing the full detector length";
  return n;
}

size_t pma::PMAlgCosmicTagger::fullWidthCrossing ( pma::TrkCandidateColl &  tracks ) const

private

Definition at line 318 of file PMAlgCosmicTagger.cxx.

{

  // Just use the first tpc to get the width dimension (instead of assuming x).
  auto const* geom = lar::providerFrom<geo::Geometry>();
  TVector3 dir = geom->TPC(0, 0).WidthDir();

  size_t n = fullCrossingTagger(tracks, ConvertDirToInt(dir));

  mf::LogInfo("pma::PMAlgCosmicTagger")
    << " - Tagged " << n << " tracks crossing the full detector width";
  return n;
}

void pma::PMAlgCosmicTagger::GetDimensions ( )

private

Definition at line 411 of file PMAlgCosmicTagger.cxx.

{

  // Need to find the minimum and maximum height values from the geometry.
  double minX = 1.e6;
  double maxX = -1.e6;
  double minY = 1.e6;
  double maxY = -1.e6;
  double minZ = 1.e6;
  double maxZ = -1.e6;

  auto const* geom = lar::providerFrom<geo::Geometry>();

  // Since we can stack TPCs, we can't just use geom::TPCGeom::Height()
  for (geo::TPCID const& tID : geom->IterateTPCIDs()) {
    geo::TPCGeo const& TPC = geom->TPC(tID);

    // We need to make sure we only include the real TPCs
    // We have dummy TPCs in the protoDUNE and DUNE geometries
    // The dummy ones have a drift distance of only ~13 cm.
    if (TPC.DriftDistance() < 25.0) { continue; }

    // get center in world coordinates
    double origin[3] = {0.};
    double center[3] = {0.};
    TPC.LocalToWorld(origin, center);
    double tpcDim[3] = {TPC.HalfWidth(), TPC.HalfHeight(), 0.5 * TPC.Length()};

    if (center[0] - tpcDim[0] < minX) minX = center[0] - tpcDim[0];
    if (center[0] + tpcDim[0] > maxX) maxX = center[0] + tpcDim[0];
    if (center[1] - tpcDim[1] < minY) minY = center[1] - tpcDim[1];
    if (center[1] + tpcDim[1] > maxY) maxY = center[1] + tpcDim[1];
    if (center[2] - tpcDim[2] < minZ) minZ = center[2] - tpcDim[2];
    if (center[2] + tpcDim[2] > maxZ) maxZ = center[2] + tpcDim[2];
  } // for all TPC

  fDimensionsMin.clear();
  fDimensionsMax.clear();
  fDimensionsMin.push_back(minX);
  fDimensionsMin.push_back(minY);
  fDimensionsMin.push_back(minZ);
  fDimensionsMax.push_back(maxX);
  fDimensionsMax.push_back(maxY);
  fDimensionsMax.push_back(maxZ);
}

bool pma::PMAlgCosmicTagger::isFrontBackVertex ( const TVector3 &  pos, double  tolerance, short int  dirIndx  ) const

private

Definition at line 399 of file PMAlgCosmicTagger.cxx.

{

  bool front = (fabs(pos[dirIndx] - fDimensionsMin[dirIndx]) < tolerance);
  bool back = (fabs(pos[dirIndx] - fDimensionsMax[dirIndx]) < tolerance);

  return front || back;
}

bool pma::PMAlgCosmicTagger::isTopVertex ( const TVector3 &  pos, double  tolerance, short int  dirIndx  ) const

private

Definition at line 392 of file PMAlgCosmicTagger.cxx.

{

  return (fabs(pos[dirIndx] - fDimensionsMax[dirIndx]) < tolerance);
}

size_t pma::PMAlgCosmicTagger::nonBeamT0Tag ( detinfo::DetectorClocksData const &  clockData, pma::TrkCandidateColl &  tracks  ) const

private

Definition at line 107 of file PMAlgCosmicTagger.cxx.

{
  size_t n = 0;

  // Search through all of the tracks
  for (auto& t : tracks.tracks()) {

    // Non zero T0 means we reconstructed it
    if (t.Track()->GetT0() != 0.0) {
      mf::LogInfo("pma::PMAlgCosmicTagger") << " - track with T0 = " << t.Track()->GetT0();

      if (fabs(t.Track()->GetT0() - trigger_offset(clockData)) > fNonBeamT0Margin) {
        ++n;
        t.Track()->SetTagFlag(pma::Track3D::kCosmic);
        t.Track()->SetTagFlag(pma::Track3D::kBeamIncompatible);
      }
    }
  }

  mf::LogInfo("pma::PMAlgCosmicTagger") << " - Tagged " << n << " non-beam T0 tracks.";
  return n;
}

size_t pma::PMAlgCosmicTagger::outOfDriftWindow ( pma::TrkCandidateColl &  tracks ) const

private

Definition at line 46 of file PMAlgCosmicTagger.cxx.

{
  mf::LogInfo("pma::PMAlgCosmicTagger") << "   - tag tracks out of 1 drift window;";
  size_t n = 0;

  auto const* geom = lar::providerFrom<geo::Geometry>();

  for (auto& t : tracks.tracks()) {

    pma::Track3D& trk = *(t.Track());

    double min, max, p;
    bool node_out_of_drift_min = false;
    bool node_out_of_drift_max = false;
    for (size_t nidx = 0; nidx < trk.Nodes().size(); ++nidx) {
      auto const& node = *(trk.Nodes()[nidx]);
      auto const& tpcGeo = geom->TPC(node.TPC(), node.Cryo());
      // DetectDriftDirection returns a short int, but switch requires an int
      int driftDir = abs(tpcGeo.DetectDriftDirection());
      p = node.Point3D()[driftDir - 1];
      switch (driftDir) {
      case 1:
        min = tpcGeo.MinX();
        max = tpcGeo.MaxX();
        break;
      case 2:
        min = tpcGeo.MinY();
        max = tpcGeo.MaxY();
        break;
      case 3:
        min = tpcGeo.MinZ();
        max = tpcGeo.MaxZ();
        break;
      default:
        throw cet::exception("PMAlgCosmicTagger")
          << "Drift direction unknown: " << driftDir << std::endl;
      }
      if (p < min - fOutOfDriftMargin) { node_out_of_drift_min = true; }
      if (p > max + fOutOfDriftMargin) { node_out_of_drift_max = true; }
    }

    if (node_out_of_drift_min && node_out_of_drift_max) {
      trk.SetTagFlag(pma::Track3D::kCosmic);
      trk.SetTagFlag(pma::Track3D::kOutsideDrift_Complete);
      ++n;
    }
    else if (node_out_of_drift_min || node_out_of_drift_max) {
      trk.SetTagFlag(pma::Track3D::kCosmic);
      trk.SetTagFlag(pma::Track3D::kOutsideDrift_Partial);
      ++n;
    }
  }

  mf::LogInfo("pma::PMAlgCosmicTagger") << " - Tagged " << n << " tracks out of 1 drift window.";

  return n;
}

void pma::PMAlgCosmicTagger::tag	(	detinfo::DetectorClocksData const &	clockData,
		pma::TrkCandidateColl &	tracks
	)

Definition at line 23 of file PMAlgCosmicTagger.cxx.

{
  // Get the detector dimensions
  GetDimensions();

  mf::LogInfo("pma::PMAlgCosmicTagger")
    << "Passed " << tracks.size() << " tracks for tagging cosmics.";

  size_t n = 0;

  if (fTagOutOfDriftTracks) n += outOfDriftWindow(tracks);
  if (fTagFullHeightTracks) n += fullHeightCrossing(tracks);
  if (fTagFullWidthTracks) n += fullWidthCrossing(tracks);
  if (fTagFullLengthTracks) n += fullLengthCrossing(tracks);
  if (fTagNonBeamT0Tracks) n += nonBeamT0Tag(clockData, tracks);
  if (fTagTopFrontBack) n += tagTopFrontBack(tracks);
  if (fTagApparentStopper) n += tagApparentStopper(tracks);

  mf::LogInfo("pma::PMAlgCosmicTagger") << "...tagged " << n << " cosmic-like tracks.";
}

bool pma::PMAlgCosmicTagger::tagAny ( ) const

inline

Definition at line 102 of file PMAlgCosmicTagger.h.

  {
    return (fTagOutOfDriftTracks || fTagFullHeightTracks || fTagFullWidthTracks ||
            fTagFullLengthTracks || fTagNonBeamT0Tracks || fTagApparentStopper || fTagTopFrontBack);
  }

size_t pma::PMAlgCosmicTagger::tagApparentStopper ( pma::TrkCandidateColl &  tracks ) const

private

Definition at line 176 of file PMAlgCosmicTagger.cxx.

{

  size_t n = 0;

  // Tracks entering from the top of the detector that stop in the fiducial volume
  // are very likely to be cosmics that leave through the APA, but have their
  // drift coordinate incorrectly set due to lack of T0
  auto const* geom = lar::providerFrom<geo::Geometry>();

  short int hIdx = ConvertDirToInt(geom->TPC(0, 0).HeightDir());

  // Loop over the tracks
  for (auto& t : tracks.tracks()) {

    // Get the first and last positions from the track.
    auto const& node0 = *(t.Track()->Nodes()[0]);
    auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);

    // Check which end is the vertex (assume the largest height)
    TVector3 vtx =
      (node0.Point3D()[hIdx] > node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
    TVector3 end =
      (node0.Point3D()[hIdx] <= node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();

    if (fabs(vtx[hIdx] - fDimensionsMax[hIdx]) < fApparentStopperMargin) {
      // Check the other element to see if it ends away from the bottom of the detector
      if (fabs(end[hIdx] - fDimensionsMin[hIdx]) > 5. * fApparentStopperMargin) {

        // We now need to loop over all of the tracks to see if any start within fStopperBuffer of our end point.
        bool foundTrack = false;
        for (auto const& tt : tracks.tracks()) {
          // Don't match with itself!
          if ((&tt) == (&t)) continue;

          // Compare this track with our main track
          TVector3 trkVtx = (tt.Track()->Nodes()[0])->Point3D();
          TVector3 trkEnd = (tt.Track()->Nodes()[tt.Track()->Nodes().size() - 1])->Point3D();

          if ((end - trkVtx).Mag() < fStopperBuffer || (end - trkEnd).Mag() < fStopperBuffer) {
            foundTrack = true;
            break;
          }
        }
        if (foundTrack) {
          // This isn't really a stopping particle, so move on
          continue;
        }

        // If we don't mind about tagging all stopping particles then this satisfies our requirements
        if (!fVetoActualStopper) {
          ++n;
          t.Track()->SetTagFlag(pma::Track3D::kCosmic);
          t.Track()->SetTagFlag(pma::Track3D::kGeometry_Y);
          continue;
        }

        // If we want to actually ignore the stopping particles, use de/dx...
        // Store the number of sigma from the mean for the final dedx point in each view
        std::vector<float> nSigmaPerView;

        // Loop over the views
        for (auto const view : geom->Views()) {

          // Get the dedx for this track and view
          std::map<size_t, std::vector<double>> dedx;
          t.Track()->GetRawdEdxSequence(dedx, view);

          std::vector<double> trk_dedx;

          for (int h = t.Track()->NextHit(-1, view); h != -1; h = t.Track()->NextHit(h, view)) {
            // If this is the last hit then this method won't work
            if (h > t.Track()->PrevHit(t.Track()->size(), view)) break;
            // Make sure we have a reasonable value
            if (dedx[h][5] / dedx[h][6] <= 0 || dedx[h][5] / dedx[h][6] > 1e6) continue;
            trk_dedx.push_back(dedx[h][5] / dedx[h][6]);
          }

          if (trk_dedx.size() == 0) {
            mf::LogInfo("pma::PMAlgCosmicTagger")
              << "View " << view << " has no hits." << std::endl;
            continue;
          }

          double sum = std::accumulate(std::begin(trk_dedx), std::end(trk_dedx), 0.0);
          double mean = sum / static_cast<double>(trk_dedx.size());
          double accum = 0.0;
          std::for_each(std::begin(trk_dedx), std::end(trk_dedx), [&](const double d) {
            accum += (d - mean) * (d - mean);
          });
          double stdev = sqrt(accum / static_cast<double>(trk_dedx.size() - 1));

          mf::LogInfo("pma::PMAlgCosmicTagger")
            << " View " << view << " has average dedx " << mean << " +/- " << stdev
            << " and final dedx " << trk_dedx[trk_dedx.size() - 1] << std::endl;

          nSigmaPerView.push_back(fabs((trk_dedx[trk_dedx.size() - 1] - mean) / stdev));
        }

        bool notStopper = true;
        short unsigned int n2Sigma = 0;
        short unsigned int n3Sigma = 0;
        for (auto const nSigma : nSigmaPerView) {
          if (nSigma >= 2.0) ++n2Sigma;
          if (nSigma >= 3.0) ++n3Sigma;
        }

        if (n3Sigma > 0) notStopper = false;
        if (n2Sigma == nSigmaPerView.size()) notStopper = false;

        if (notStopper) {
          ++n;
          t.Track()->SetTagFlag(pma::Track3D::kCosmic);
          t.Track()->SetTagFlag(pma::Track3D::kGeometry_Y);
        }
      } // Check on bottom position
    }   // Check on top position
  }     // End loop over tracks

  mf::LogInfo("pma::PMAlgCosmicTagger")
    << " - Tagged " << n << " tracks stopping in the detector after starting at the top."
    << std::endl;

  return n;
}

size_t pma::PMAlgCosmicTagger::tagTopFrontBack ( pma::TrkCandidateColl &  tracks ) const

private

Definition at line 132 of file PMAlgCosmicTagger.cxx.

{

  size_t n = 0;

  auto const* geom = lar::providerFrom<geo::Geometry>();

  short int hIdx = ConvertDirToInt(geom->TPC(0, 0).HeightDir());
  short int lIdx = ConvertDirToInt(geom->TPC(0, 0).LengthDir());

  // Loop over the tracks
  for (auto& t : tracks.tracks()) {

    // Get the first and last positions from the track.
    auto const& node0 = *(t.Track()->Nodes()[0]);
    auto const& node1 = *(t.Track()->Nodes()[t.Track()->Nodes().size() - 1]);

    // Check which end is the vertex (assume the largest height)
    TVector3 vtx =
      (node0.Point3D()[hIdx] > node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();
    TVector3 end =
      (node0.Point3D()[hIdx] <= node1.Point3D()[hIdx]) ? node0.Point3D() : node1.Point3D();

    // Check we have a track starting at the top of the detector
    bool top = isTopVertex(vtx, fTopFrontBackMargin, hIdx);

    // Check the track ends at the front or back of the detector
    bool frontBack = isFrontBackVertex(end, fTopFrontBackMargin, lIdx);

    // Check we path both criteria but without letting either the start or end of the track fulfill both
    if (top && frontBack) {
      ++n;
      t.Track()->SetTagFlag(pma::Track3D::kCosmic);
      t.Track()->SetTagFlag(pma::Track3D::kGeometry_YZ);
    }
  }

  mf::LogInfo("pma::PMAlgCosmicTagger")
    << " - Tagged " << n << " tracks crossing from top to front/back." << std::endl;

  return n;
}

Member Data Documentation

double pma::PMAlgCosmicTagger::fApparentStopperMargin

private

Definition at line 146 of file PMAlgCosmicTagger.h.

std::vector<double> pma::PMAlgCosmicTagger::fDimensionsMax

private

Definition at line 155 of file PMAlgCosmicTagger.h.

std::vector<double> pma::PMAlgCosmicTagger::fDimensionsMin

private

Definition at line 154 of file PMAlgCosmicTagger.h.

double pma::PMAlgCosmicTagger::fFullCrossingMargin

private

Definition at line 136 of file PMAlgCosmicTagger.h.

double pma::PMAlgCosmicTagger::fNonBeamT0Margin

private

Definition at line 140 of file PMAlgCosmicTagger.h.

double pma::PMAlgCosmicTagger::fOutOfDriftMargin

private

Definition at line 130 of file PMAlgCosmicTagger.h.

double pma::PMAlgCosmicTagger::fStopperBuffer

private

Definition at line 148 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fTagApparentStopper

private

Definition at line 145 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fTagFullHeightTracks

private

Definition at line 133 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fTagFullLengthTracks

private

Definition at line 135 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fTagFullWidthTracks

private

Definition at line 134 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fTagNonBeamT0Tracks

private

Definition at line 139 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fTagOutOfDriftTracks

private

Definition at line 129 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fTagTopFrontBack

private

Definition at line 142 of file PMAlgCosmicTagger.h.

double pma::PMAlgCosmicTagger::fTopFrontBackMargin

private

Definition at line 143 of file PMAlgCosmicTagger.h.

bool pma::PMAlgCosmicTagger::fVetoActualStopper

private

Definition at line 147 of file PMAlgCosmicTagger.h.

---

The documentation for this class was generated from the following files:

• PMAlgCosmicTagger.h
• PMAlgCosmicTagger.cxx