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 "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// LArSoft includes
#include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom<>()
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "larcore/Geometry/Geometry.h"

// ROOT
#include <TTree.h>
#include <TH2F.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TCutG.h>
#include <TString.h>
#include <TMarker.h>
#include <TColor.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TVirtualPad.h>
#include <TLatex.h>
#include <TGraph.h>
#include <TF1.h>
#include <TLine.h>
#include <TPrincipal.h>
#include <TMath.h>
#include <TVector.h>
#include <TVectorD.h>
#include <TVector2.h>

// c++
#include <string>
#include <vector>
#include <map>
#include <sstream>


namespace cluster {
  class BlurredClusteringAlg;
}

class cluster::BlurredClusteringAlg {
public:

  BlurredClusteringAlg(fhicl::ParameterSet const& pset);
  virtual ~BlurredClusteringAlg();

  void reconfigure(fhicl::ParameterSet const&p);

  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);

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

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

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

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

  unsigned int GetMinSize() { return fMinSize; }

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

  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);

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

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

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

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

  void FindBlurringParameters(int& blurwire, int& blurtick, int& sigmawire, int& sigmatick);

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

  void MakeKernels();

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

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

  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 fMinMergeClusterSize;      // minimum size of a cluster to consider merging it to another
  unsigned int fNeighboursThreshold;      // min. number of neighbors to add to cluster
  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;
  std::string fDebugPDFName;

  // art service handles
  art::ServiceHandle<geo::Geometry> fGeom;
  detinfo::DetectorProperties const* fDetProp;
  lariov::ChannelStatusProvider const& fChanStatus = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();

};

#endif