HoughVertex_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief  Find event vertex using a Hough transform
/// \author Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "fhiclcpp/ParameterSet.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/Vertex.h"
#include "Utilities/AssociationUtil.h"

#include "HoughVertex/HoughVertexAlg.h"

#include "TH2.h"

#include <string>

namespace hv
{
  /// Find event vertex using a Hough transform
  class HoughVertex: public art::EDProducer
  {
  public:
    explicit HoughVertex(const fhicl::ParameterSet& pset);
    virtual ~HoughVertex();

    void produce(art::Event& evt);

  protected:
    std::string fSliceLabel;
    double fErrorScale;
    double fMaxHitsUsedPerView;
    bool fObeyPreselection;

    HoughVertexAlg fAlg;

    bool fSaveMaps;
  };

  //......................................................................
  HoughVertex::HoughVertex(const fhicl::ParameterSet& pset)
    : fSliceLabel(pset.get<std::string>("SliceLabel")),
      fErrorScale(pset.get<double>("ErrorScale")),
      fMaxHitsUsedPerView(pset.get<double>("MaxHitsUsedPerView")),
      fObeyPreselection(pset.get<bool>("ObeyPreselection")),
      fAlg(pset.get<fhicl::ParameterSet>("VertexAlgPSet")),
      fSaveMaps(pset.get<bool>("SaveMaps"))
  {
    produces<std::vector<rb::Vertex>>();
    // Vertices are associated with the slice they came from
    produces<art::Assns<rb::Vertex, rb::Cluster>>();

    if(fSaveMaps) produces<std::vector<TH2F>>();
  }

  //......................................................................
  HoughVertex::~HoughVertex()
  {
  }

  //......................................................................
  void HoughVertex::produce(art::Event& evt)
  {
    std::unique_ptr<std::vector<rb::Vertex>> vtxcol(new std::vector<rb::Vertex>);
    std::unique_ptr<std::vector<TH2F>> histcol(new std::vector<TH2F>);
    // Associate vertices with the slice they came from
    std::unique_ptr<art::Assns<rb::Vertex, rb::Cluster>> assns(new art::Assns<rb::Vertex, rb::Cluster>);

    art::ServiceHandle<geo::Geometry> geom;

    art::Handle<std::vector<rb::Cluster>> slices;
    evt.getByLabel(fSliceLabel, slices);

    for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
      if(fObeyPreselection && rb::IsFiltered(evt, slices, sliceIdx)) continue;

      const rb::Cluster& slice = (*slices)[sliceIdx];
      if(slice.IsNoise()) continue;

      std::vector<Vector2D> xhits, yhits;

      for(unsigned int hitIdx = 0; hitIdx < slice.NCell(); ++hitIdx){
        const rb::CellHit& chit = *slice.Cell(hitIdx);

        TVector3 xyz;
        const geo::CellGeo* geocell = geom->Plane(chit.Plane())->Cell(chit.Cell());
        geocell->GetCenter(xyz);
        const float dz = (2*geocell->HalfD())/sqrt(12) * fErrorScale;
        const float dv = (2*geocell->HalfW())/sqrt(12) * fErrorScale;

        if(chit.View() == geo::kX){
          xhits.push_back(Vector2D(xyz.Z(), xyz.X(), dz, dv));
        }
        else{
          yhits.push_back(Vector2D(xyz.Z(), xyz.Y(), dz, dv));
        }
      } // end for hitIdx

      // Vertexing time goes quadratically as the number of hits. Only use the
      // N most-upstream hits in each view.
      std::sort(xhits.begin(), xhits.end());
      std::sort(yhits.begin(), yhits.end());
      while(xhits.size() > fMaxHitsUsedPerView) xhits.pop_back();
      while(yhits.size() > fMaxHitsUsedPerView) yhits.pop_back();

      // Do all the actual work
      TVector3 vtx;
      if(fSaveMaps){
        TH2F *xmap, *ymap;

        vtx = fAlg.FindVertex(xhits, yhits, &xmap, &ymap);

        xmap->SetName(TString::Format("xmap_slice%d", sliceIdx));
        ymap->SetName(TString::Format("ymap_slice%d", sliceIdx));

        xmap->SetTitle(TString::Format("X view Hough map, slice %d", sliceIdx));
        ymap->SetTitle(TString::Format("Y view Hough map, slice %d", sliceIdx));

        histcol->push_back(*xmap);
        histcol->push_back(*ymap);

        delete xmap;
        delete ymap;
      }
      else{
        vtx = fAlg.FindVertex(xhits, yhits);
      }

      vtxcol->push_back(rb::Vertex(vtx, slice.MeanTNS()));
      // Associate most recently added vertex with the slice
      util::CreateAssn(*this, evt, *vtxcol,
                       art::Ptr<rb::Cluster>(slices, sliceIdx), *assns);
    } // end for sliceIdx

    evt.put(std::move(vtxcol));
    if(fSaveMaps) evt.put(std::move(histcol));
    evt.put(std::move(assns));
  }

  DEFINE_ART_MODULE(HoughVertex)

} // namespace
//////////////////////////////////////////////////////////////////////////