quad Namespace Reference

Classes
class	EvalVtx
class	HeatMap
struct	Line2D
struct	Pt2D
class	QuadVtx

Functions
recob::Vertex	GetFirstVertex (const std::string &label, const art::Event &evt)
recob::Vertex	GetVtxByAssns (const std::string &label, const art::Event &evt)
bool	IntersectsCircle (float m, float c, float z0, float x0, float R, float &z1, float &z2)
void	LinesFromPoints (const std::vector< Pt2D > &pts, std::vector< Line2D > &lines, float z0=0, float x0=0, float R=-1)
bool	CloseAngles (float ma, float mb)
void	MapFromLines (const std::vector< Line2D > &lines, HeatMap &hm)
TVector3	FindPeak3D (const std::vector< HeatMap > &hs, const std::vector< TVector3 > &dirs) noexcept
void	GetPts2D (const detinfo::DetectorPropertiesData &detProp, const std::vector< recob::Hit > &hits, std::vector< std::vector< Pt2D >> &pts, std::vector< TVector3 > &dirs, const geo::GeometryCore *geom)

Function Documentation

bool quad::CloseAngles ( float  ma, float  mb  )

inline

Definition at line 188 of file QuadVtx_module.cc.

  {
    const float cosCrit = cos(10 * M_PI / 180.);
    const float dot = 1 + ma * mb; // (1, ma)*(1, mb)
    return cet::square(dot) > (1 + cet::square(ma)) * (1 + cet::square(mb)) * cet::square(cosCrit);
  }

TVector3 quad::FindPeak3D ( const std::vector< HeatMap > &  hs, const std::vector< TVector3 > &  dirs  )

noexcept

Definition at line 259 of file QuadVtx_module.cc.

  {
    assert(hs.size() == 3);
    assert(dirs.size() == 3);

    const int Nx = hs[0].Nx;

    TMatrixD M(2, 2);
    M(0, 0) = dirs[0].Y();
    M(0, 1) = dirs[0].Z();
    M(1, 0) = dirs[1].Y();
    M(1, 1) = dirs[1].Z();

    // Singular, and stupid setup of exceptions means we can't test any other way
    if (M(0, 0) * M(1, 1) - M(1, 0) * M(0, 1) == 0) return TVector3(0, 0, 0);

    M.Invert();

    float bestscore = -1;
    TVector3 bestr;

    // Accumulate some statistics up front that will enable us to optimize
    std::vector<float> colMax[3];
    for (int view = 0; view < 3; ++view) {
      colMax[view].resize(hs[view].Nz);
      for (int iz = 0; iz < hs[view].Nz; ++iz) {
        colMax[view][iz] = *std::max_element(&hs[view].map[Nx * iz], &hs[view].map[Nx * (iz + 1)]);
      }
    }

    for (int iz = 0; iz < hs[0].Nz; ++iz) {
      const float z = hs[0].ZBinCenter(iz);
      const float bonus = 1; // works badly... exp((hs[0].maxz-z)/1000.);

      for (int iu = 0; iu < hs[1].Nz; ++iu) {
        const float u = hs[1].ZBinCenter(iu);
        // r.Dot(d0) = z && r.Dot(d1) = u
        TVectorD p(2);
        p(0) = z;
        p(1) = u;
        const TVectorD r = M * p;
        const float v = r[0] * dirs[2].Y() + r[1] * dirs[2].Z();
        const int iv = hs[2].ZToBin(v);
        if (iv < 0 || iv >= hs[2].Nz) continue;
        const double y = r(0);

        // Even if the maxes were all at the same x we couldn't beat the record
        if (colMax[0][iz] + colMax[1][iu] + colMax[2][iv] < bestscore) continue;

        // Attempt to micro-optimize the dx loop below
        const float* __restrict__ h0 = &hs[0].map[Nx * iz];
        const float* __restrict__ h1 = &hs[1].map[Nx * iu];
        const float* __restrict__ h2 = &hs[2].map[Nx * iv];

        int bestix = -1;
        for (int ix = 1; ix < Nx - 1; ++ix) {
          const float score = bonus * (h0[ix] + h1[ix] + h2[ix]);

          if (score > bestscore) {
            bestscore = score;
            bestix = ix;
          }
        } // end for dx

        if (bestix != -1) { bestr = TVector3(hs[0].XBinCenter(bestix), y, z); }
      } // end for u
    }   // end for z

    return bestr;
  }

recob::Vertex quad::GetFirstVertex	(	const std::string &	label,
		const art::Event &	evt
	)

Definition at line 70 of file EvalVtx_module.cc.

{
  const auto& vtxs = *evt.getValidHandle<std::vector<recob::Vertex>>(label);

  if(vtxs.empty()) return recob::Vertex();

  return vtxs[0];
}

void quad::GetPts2D	(	const detinfo::DetectorPropertiesData &	detProp,
		const std::vector< recob::Hit > &	hits,
		std::vector< std::vector< Pt2D >> &	pts,
		std::vector< TVector3 > &	dirs,
		const geo::GeometryCore *	geom
	)

Definition at line 332 of file QuadVtx_module.cc.

  {
    pts.resize(3); // 3 views

    TVector3 dirZ(0, 0, 1);
    TVector3 dirU, dirV;

    for (const recob::Hit& hit : hits) {
      const geo::WireID wire = hit.WireID();

      const double xpos = detProp.ConvertTicksToX(hit.PeakTime(), wire);

      const TVector3 r0 = geom->WireEndPoints(wire).start();
      const TVector3 r1 = geom->WireEndPoints(wire).end();

      const double energy = hit.Integral();

      if (geom->View(hit.Channel()) == geo::kZ) {
        pts[0].emplace_back(xpos, r0.z(), 0, energy);
        continue;
      }

      // Compute the direction perpendicular to the wires
      TVector3 perp = (r1 - r0).Unit();
      perp = TVector3(0, -perp.z(), perp.y());
      // We want to ultimately have a positive z component in "perp"
      if (perp.z() < 0) perp *= -1;

      // TODO check we never get a 4th view a-la the bug we had in the 3D version

      // The "U" direction is the first one we see
      if (dirU.Mag2() == 0) { dirU = perp; }
      else if (dirV.Mag2() == 0 && fabs(dirU.Dot(perp)) < 0.99) {
        // If we still need a "V" and this direction differs from "U"
        dirV = perp;
      }

      // Hits belong to whichever view their perpendicular vector aligns with
      if (fabs(dirU.Dot(perp)) > 0.99) { pts[1].emplace_back(xpos, r0.Dot(dirU), 1, energy); }
      else {
        pts[2].emplace_back(xpos, r0.Dot(dirV), 2, energy);
      }
    } // end for hits

    dirs = {dirZ, dirU, dirV};

    std::default_random_engine gen;

    // In case we need to sub-sample they should be shuffled
    for (int view = 0; view < 3; ++view) {
      std::shuffle(pts[view].begin(), pts[view].end(), gen);
    }
  }

recob::Vertex quad::GetVtxByAssns	(	const std::string &	label,
		const art::Event &	evt
	)

Definition at line 80 of file EvalVtx_module.cc.

{
  art::Handle<std::vector<recob::Vertex>> vtxs;
  evt.getByLabel(label, vtxs);

  if(vtxs->empty()) return recob::Vertex();

  art::Handle<std::vector<recob::PFParticle>> parts;
  evt.getByLabel(label, parts);

  art::FindManyP<recob::Vertex> fm(parts, evt, label);

  for(unsigned int i = 0; i < parts->size(); ++i){
    const int pdg = abs((*parts)[i].PdgCode());
    if(pdg == 12 || pdg == 14 || pdg == 16){
      const std::vector<art::Ptr<recob::Vertex>>& vtxs = fm.at(i);
      if(vtxs.size() == 1) return *vtxs[0];

      if(vtxs.empty()){
        std::cout << "Warning: vertex list empty!" << std::endl;
      }
      else{
        std::cout << "Warning: " << vtxs.size() << " vertices for daughter?" << std::endl;
      }
    }
  }

  return recob::Vertex();
}

bool quad::IntersectsCircle	(	float	m,
		float	c,
		float	z0,
		float	x0,
		float	R,
		float &	z1,
		float &	z2
	)

Definition at line 115 of file QuadVtx_module.cc.

  {
    // Change to the frame where (z0, x0) = (0, 0)
    c += m * z0 - x0;

    // z^2 + (m*z+c)^2 = R^2
    const float A = 1 + cet::square(m);
    const float B = 2 * m * c;
    const float C = cet::square(c) - cet::square(R);

    double desc = cet::square(B) - 4 * A * C;

    if (desc < 0) return false;

    desc = sqrt(desc);

    z1 = (-B - desc) / (2 * A);
    z2 = (-B + desc) / (2 * A);

    // Back to the original frame
    z1 += z0;
    z2 += z0;

    return true;
  }

void quad::LinesFromPoints	(	const std::vector< Pt2D > &	pts,
		std::vector< Line2D > &	lines,
		float	z0 = 0,
		float	x0 = 0,
		float	R = -1
	)

Definition at line 143 of file QuadVtx_module.cc.

  {
    constexpr size_t kMaxLines = 10 * 1000 * 1000; // This is 150MB of lines...

    const size_t product = (pts.size() * (pts.size() - 1)) / 2;
    const int stride = product / kMaxLines + 1;

    lines.reserve(std::min(product, kMaxLines));

    for (int offset = 0; offset < stride; ++offset) {
      for (unsigned int i = 0; i < pts.size(); ++i) {
        for (unsigned int j = i + offset + 1; j < pts.size(); j += stride) {
          const Line2D l(pts[i], pts[j]);

          if (isinf(l.m) || isnan(l.m) || isinf(l.c) || isnan(l.c)) continue;

          if (R > 0) {
            float z1, z2;
            if (!IntersectsCircle(l.m, l.c, z0, x0, 2.5, z1, z2)) continue;
            if (l.minz < z1 && l.minz < z2 && l.maxz > z1 && l.maxz > z2) continue;
          }

          lines.push_back(std::move(l));
          if (lines.size() == kMaxLines) goto end; // break out of 3 loops
        }
      }
    }

  end:

    lines.shrink_to_fit();

    mf::LogInfo() << "Made " << lines.size() << " lines using stride " << stride
                  << " to fit under cap of " << kMaxLines << std::endl;

    // Lines are required to be sorted by gradient for a later optimization
    std::sort(lines.begin(), lines.end());
  }

void quad::MapFromLines	(	const std::vector< Line2D > &	lines,
		HeatMap &	hm
	)

Definition at line 197 of file QuadVtx_module.cc.

  {
    // This maximum is driven by runtime
    constexpr size_t kMaxPts = 10 * 1000 * 1000;

    unsigned int j0 = 0;
    unsigned int jmax = 0;

    long npts = 0;
    for (unsigned int i = 0; i + 1 < lines.size(); ++i) {
      const Line2D a = lines[i];

      j0 = std::max(j0, i + 1);
      while (j0 < lines.size() && CloseAngles(a.m, lines[j0].m))
        ++j0;
      jmax = std::max(jmax, j0);
      while (jmax < lines.size() && !CloseAngles(a.m, lines[jmax].m))
        ++jmax;

      npts += jmax - j0;
    }

    const size_t product = (lines.size() * (lines.size() - 1)) / 2;
    const int stride = npts / kMaxPts + 1;

    mf::LogInfo() << "Combining lines to points with stride " << stride << std::endl;

    mf::LogInfo() << npts << " cf " << product << " ie " << double(npts) / product << std::endl;

    j0 = 0;
    jmax = 0;

    for (unsigned int i = 0; i + 1 < lines.size(); ++i) {
      const Line2D a = lines[i];

      j0 = std::max(j0, i + 1);
      while (j0 < lines.size() && CloseAngles(a.m, lines[j0].m))
        ++j0;
      jmax = std::max(jmax, j0);
      while (jmax < lines.size() && !CloseAngles(a.m, lines[jmax].m))
        ++jmax;

      for (unsigned int j = j0; j < jmax; j += stride) {
        const Line2D& b = lines[j];

        // x = mA * z + cA = mB * z + cB
        const float z = (b.c - a.c) / (a.m - b.m);
        const float x = a.m * z + a.c;

        // No solutions within a line
        if ((z < a.minz || z > a.maxz) && (z < b.minz || z > b.maxz)) {
          const int iz = hm.ZToBin(z);
          const int ix = hm.XToBin(x);
          if (iz >= 0 && iz < hm.Nz && ix >= 0 && ix < hm.Nx) { hm.map[iz * hm.Nx + ix] += stride; }
        }
      } // end for i
    }
  }