quad::QuadVtx Class Reference

Inheritance diagram for quad::QuadVtx:

Public Member Functions
	QuadVtx (const fhicl::ParameterSet &pset)
void	beginJob () override
void	produce (art::Event &evt) override

Private Member Functions
bool	FindVtx (const detinfo::DetectorPropertiesData &detProp, const std::vector< recob::Hit > &hits, TVector3 &vtx, int evt) const

Private Attributes
std::string	fHitLabel
bool	fSavePlots
const geo::GeometryCore *	geom

Detailed Description

Definition at line 74 of file QuadVtx_module.cc.

Constructor & Destructor Documentation

quad::QuadVtx::QuadVtx ( const fhicl::ParameterSet &  pset )

explicit

Definition at line 97 of file QuadVtx_module.cc.

    : EDProducer(pset)
    , fHitLabel(pset.get<std::string>("HitLabel"))
    , fSavePlots(pset.get<bool>("SavePlots"))
  {
    produces<std::vector<recob::Vertex>>();
  }

Member Function Documentation

void quad::QuadVtx::beginJob ( )

override

Definition at line 107 of file QuadVtx_module.cc.

  {
    geom = art::ServiceHandle<const geo::Geometry>()->provider();
  }

bool quad::QuadVtx::FindVtx ( const detinfo::DetectorPropertiesData &  detProp, const std::vector< recob::Hit > &  hits, TVector3 &  vtx, int  evt  ) const

private

Definition at line 392 of file QuadVtx_module.cc.

  {
    if (hits.empty()) return false;

    std::vector<std::vector<Pt2D>> pts;
    std::vector<TVector3> dirs;

    GetPts2D(detProp, hits, pts, dirs, geom);

    double minx = +1e9;
    double maxx = -1e9;
    double minz[3] = {+1e9, +1e9, +1e9};
    double maxz[3] = {-1e9, -1e9, -1e9};
    for (int view = 0; view < 3; ++view) {
      for (const Pt2D& p : pts[view]) {
        minx = std::min(minx, p.x);
        maxx = std::max(maxx, p.x);
        minz[view] = std::min(minz[view], p.z);
        maxz[view] = std::max(maxz[view], p.z);
      }
    }

    // Add some padding
    for (int view = 0; view < 3; ++view) {
      minz[view] -= 100;
      maxz[view] += 100;
    }
    minx -= 20;
    maxx += 20;

    // Don't allow the vertex further downstream in z (view 0) than 25% of the
    // hits.
    std::vector<float> zs;
    zs.reserve(pts[0].size());
    for (const Pt2D& p : pts[0])
      zs.push_back(p.z);
    auto mid = zs.begin() + zs.size() / 4;
    if (mid != zs.end()) {
      std::nth_element(zs.begin(), mid, zs.end());
      maxz[0] = *mid;
    }

    std::vector<HeatMap> hms;
    hms.reserve(3);
    for (int view = 0; view < 3; ++view) {
      if (pts[view].empty()) return false;

      std::vector<Line2D> lines;
      LinesFromPoints(pts[view], lines);

      if (lines.empty()) return false;

      // Approximately cm bins
      hms.emplace_back(maxz[view] - minz[view], minz[view], maxz[view], maxx - minx, minx, maxx);
      MapFromLines(lines, hms.back());
    } // end for view

    vtx = FindPeak3D(hms, dirs);

    std::vector<HeatMap> hms_zoom;
    hms_zoom.reserve(3);
    for (int view = 0; view < 3; ++view) {
      const double x0 = vtx.X();
      const double z0 = vtx.Dot(dirs[view]);

      std::vector<Line2D> lines;
      LinesFromPoints(pts[view], lines, z0, x0, 2.5);

      if (lines.empty()) return false; // How does this happen??

      // mm granularity
      hms_zoom.emplace_back(50, z0 - 2.5, z0 + 2.5, 50, x0 - 2.5, x0 + 2.5);

      MapFromLines(lines, hms_zoom.back());
    }

    vtx = FindPeak3D(hms_zoom, dirs);

    if (fSavePlots) {
      art::TFileDirectory evt_dir =
        art::ServiceHandle<art::TFileService>()->mkdir(TString::Format("evt%d", evt).Data());

      for (int view = 0; view < 3; ++view) {
        art::TFileDirectory view_dir = evt_dir.mkdir(TString::Format("view%d", view).Data());

        TGraph* gpts = view_dir.makeAndRegister<TGraph>("hits", "");
        for (const Pt2D& p : pts[view])
          gpts->SetPoint(gpts->GetN(), p.z, p.x);

        view_dir.makeAndRegister<TH2F>("hmap", "", *hms[view].AsTH2());

        view_dir.makeAndRegister<TH2F>("hmap_zoom", "", *hms_zoom[view].AsTH2());

        const double x = vtx.X();
        const double z = vtx.Dot(dirs[view]);
        view_dir.makeAndRegister<TGraph>("vtx3d", "", 1, &z, &x);
      } // end for view
    }   // end if saving plots

    return true;
  }

void quad::QuadVtx::produce ( art::Event &  evt )

override

Definition at line 499 of file QuadVtx_module.cc.

  {
    auto vtxcol = std::make_unique<std::vector<recob::Vertex>>();

    art::Handle<std::vector<recob::Hit>> hits;
    evt.getByLabel(fHitLabel, hits);

    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt);
    TVector3 vtx;
    if (FindVtx(detProp, *hits, vtx, evt.event())) {
      vtxcol->emplace_back(
        recob::Vertex::Point_t(vtx.X(), vtx.Y(), vtx.Z()), recob::Vertex::SMatrixSym33(), 0, 0);
    }

    evt.put(std::move(vtxcol));
  }

Member Data Documentation

std::string quad::QuadVtx::fHitLabel

private

Definition at line 87 of file QuadVtx_module.cc.

bool quad::QuadVtx::fSavePlots

private

Definition at line 89 of file QuadVtx_module.cc.

const geo::GeometryCore* quad::QuadVtx::geom

private

Definition at line 91 of file QuadVtx_module.cc.

---

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

• QuadVtx_module.cc