BlurredClusteringAlg.h

Go to the documentation of this file.

// Implementation of the Blurred Clustering algorithm
//
// Converts a hit map into a 2D image of the hits before convoling
// with a Gaussian function to introduce a weighted blurring.
// Clustering proceeds on this blurred image to create more
// complete clusters.
//
// M Wallbank (m.wallbank@sheffield.ac.uk), May 2015

#ifndef BlurredClustering_h
#define BlurredClustering_h

// Framework includes
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "fhiclcpp/fwd.h"

// LArSoft includes
#include "larcore/Geometry/Geometry.h"
#include "larcore/Geometry/WireReadout.h"
#include "lardataobj/RecoBase/Hit.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
namespace detinfo {
  class DetectorProperties;
}
namespace lariov {
  class ChannelStatusProvider;
}
namespace geo {
  struct WireID;
}

// ROOT
#include "TString.h"
class TCanvas;
class TH2F;

// c++
#include <array>
#include <string>
#include <vector>

namespace cluster {
  class BlurredClusteringAlg;
}

class cluster::BlurredClusteringAlg {
public:
  BlurredClusteringAlg(fhicl::ParameterSet const& pset);
  ~BlurredClusteringAlg();

  void CreateDebugPDF(int run, int subrun, int event);

  void ConvertBinsToClusters(std::vector<std::vector<double>> const& image,
                             std::vector<std::vector<int>> const& allClusterBins,
                             std::vector<art::PtrVector<recob::Hit>>& clusters) const;

  std::vector<std::vector<double>> ConvertRecobHitsToVector(
    std::vector<art::Ptr<recob::Hit>> const& hits,
    int readoutWindowSize);

  int FindClusters(std::vector<std::vector<double>> const& image,
                   std::vector<std::vector<int>>& allcluster) const;

  int GlobalWire(geo::WireID const& wireID) const;

  std::vector<std::vector<double>> GaussianBlur(
    std::vector<std::vector<double>> const& image) const;

  unsigned int GetMinSize() const noexcept { return fMinSize; }

  TH2F* MakeHistogram(std::vector<std::vector<double>> const& image, TString name) const;

  void SaveImage(TH2F* image,
                 std::vector<art::PtrVector<recob::Hit>> const& allClusters,
                 int pad,
                 int tpc,
                 int plane);

  void SaveImage(TH2F* image, int pad, int tpc, int plane);

  void SaveImage(TH2F* image,
                 std::vector<std::vector<int>> const& allClusterBins,
                 int pad,
                 int tpc,
                 int plane);

private:
  art::PtrVector<recob::Hit> ConvertBinsToRecobHits(std::vector<std::vector<double>> const& image,
                                                    std::vector<int> const& bins) const;

  art::Ptr<recob::Hit> ConvertBinToRecobHit(std::vector<std::vector<double>> const& image,
                                            int bin) const;

  int ConvertWireTickToBin(std::vector<std::vector<double>> const& image, int xbin, int ybin) const;

  double ConvertBinToCharge(std::vector<std::vector<double>> const& image, int bin) const;

  std::pair<int, int> DeadWireCount(int wire_bin, int width) const;

  std::array<int, 4> FindBlurringParameters() const;

  double GetTimeOfBin(std::vector<std::vector<double>> const& image, int bin) const;

  std::vector<std::vector<std::vector<double>>> MakeKernels() const;

  unsigned int NumNeighbours(int nx, std::vector<bool> const& used, int bin) const;

  bool PassesTimeCut(std::vector<double> const& times, double time) const;

  bool fDebug;
  std::string fDetector;

  // Parameters used in the Blurred Clustering algorithm
  int fBlurWire;            // blur radius for Gauss kernel in the wire direction
  int fBlurTick;            // blur radius for Gauss kernel in the tick direction
  double fSigmaWire;        // sigma for Gaussian kernel in the wire direction
  double fSigmaTick;        // sigma for Gaussian kernel in the tick direction
  int fMaxTickWidthBlur;    // maximum distance to blur a hit based on its natural width in time
  int fClusterWireDistance; // how far to cluster from seed in wire direction
  int fClusterTickDistance; // how far to cluster from seed in tick direction
  unsigned int fNeighboursThreshold; // min. number of neighbors to add to cluster
  unsigned int fMinNeighbours;       // minumum number of neighbors to keep in the cluster
  unsigned int fMinSize;             // minimum size for cluster
  double fMinSeed;         // minimum seed after blurring needed before clustering proceeds
  double fTimeThreshold;   // time threshold for clustering
  double fChargeThreshold; // charge threshold for clustering

  // Blurring stuff
  int fKernelWidth, fKernelHeight;
  std::vector<std::vector<std::vector<double>>> fAllKernels;

  // Hit containers
  std::vector<std::vector<art::Ptr<recob::Hit>>> fHitMap;
  std::vector<bool> fDeadWires;

  int fLowerTick, fUpperTick;
  int fLowerWire, fUpperWire;

  // For the debug pdf
  TCanvas* fDebugCanvas{nullptr};
  std::string fDebugPDFName{};

  // art service handles
  geo::WireReadoutGeom const* fWireReadoutGeom{
    &art::ServiceHandle<geo::WireReadout const>()->Get()};
  lariov::ChannelStatusProvider const& fChanStatus{
    art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider()};
};

#endif