cluster::BlurredClusteringAlg Class Reference

#include <BlurredClusteringAlg.h>

Public Member Functions
	BlurredClusteringAlg (fhicl::ParameterSet const &pset)
	~BlurredClusteringAlg ()
void	CreateDebugPDF (int run, int subrun, int event)
	Create the PDF to save debug images. More...
void	ConvertBinsToClusters (std::vector< std::vector< double >> const &image, std::vector< std::vector< int >> const &allClusterBins, std::vector< art::PtrVector< recob::Hit >> &clusters) const
	Takes a vector of clusters (itself a vector of hits) and turns them into clusters using the initial hit selection. More...
std::vector< std::vector< double > >	ConvertRecobHitsToVector (std::vector< art::Ptr< recob::Hit >> const &hits, int readoutWindowSize)
	Takes hit map and returns a 2D vector representing wire and tick, filled with the charge. More...
int	FindClusters (std::vector< std::vector< double >> const &image, std::vector< std::vector< int >> &allcluster) const
	Find clusters in the histogram. More...
int	GlobalWire (geo::WireID const &wireID) const
	Find the global wire position. More...
std::vector< std::vector< double > >	GaussianBlur (std::vector< std::vector< double >> const &image) const
	Applies Gaussian blur to image. More...
unsigned int	GetMinSize () const noexcept
	Minimum size of cluster to save. More...
TH2F *	MakeHistogram (std::vector< std::vector< double >> const &image, TString name) const
	Converts a 2D vector in a histogram for the debug pdf. More...
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)
	Save the images for debugging. More...
void	SaveImage (TH2F *image, std::vector< std::vector< int >> const &allClusterBins, int pad, int tpc, int plane)

Private Member Functions
art::PtrVector< recob::Hit >	ConvertBinsToRecobHits (std::vector< std::vector< double >> const &image, std::vector< int > const &bins) const
	Converts a vector of bins into a hit selection - not all the hits in the bins vector are real hits. More...
art::Ptr< recob::Hit >	ConvertBinToRecobHit (std::vector< std::vector< double >> const &image, int bin) const
	Converts a bin into a recob::Hit (not all of these bins correspond to recob::Hits - some are fake hits created by the blurring) More...
int	ConvertWireTickToBin (std::vector< std::vector< double >> const &image, int xbin, int ybin) const
	Converts an xbin and a ybin to a global bin number. More...
double	ConvertBinToCharge (std::vector< std::vector< double >> const &image, int bin) const
	Returns the charge stored in the global bin value. More...
std::pair< int, int >	DeadWireCount (int wire_bin, int width) const
std::array< int, 4 >	FindBlurringParameters () const
	Dynamically find the blurring radii and Gaussian sigma in each dimension. More...
double	GetTimeOfBin (std::vector< std::vector< double >> const &image, int bin) const
	Returns the hit time of a hit in a particular bin. More...
std::vector< std::vector< std::vector< double > > >	MakeKernels () const
	Makes all the kernels which could be required given the tuned parameters. More...
unsigned int	NumNeighbours (int nx, std::vector< bool > const &used, int bin) const
	Determines the number of clustered neighbours of a hit. More...
bool	PassesTimeCut (std::vector< double > const &times, double time) const
	Determine if a hit is within a time threshold of any other hits in a cluster. More...

Private Attributes
bool	fDebug
std::string	fDetector
int	fBlurWire
int	fBlurTick
double	fSigmaWire
double	fSigmaTick
int	fMaxTickWidthBlur
int	fClusterWireDistance
int	fClusterTickDistance
unsigned int	fNeighboursThreshold
unsigned int	fMinNeighbours
unsigned int	fMinSize
double	fMinSeed
double	fTimeThreshold
double	fChargeThreshold
int	fKernelWidth
int	fKernelHeight
std::vector< std::vector< std::vector< double > > >	fAllKernels
std::vector< std::vector< art::Ptr< recob::Hit > > >	fHitMap
std::vector< bool >	fDeadWires
int	fLowerTick
int	fUpperTick
int	fLowerWire
int	fUpperWire
TCanvas *	fDebugCanvas {nullptr}
std::string	fDebugPDFName {}
art::ServiceHandle< geo::Geometry const >	fGeom
lariov::ChannelStatusProvider const &	fChanStatus

Detailed Description

Definition at line 49 of file BlurredClusteringAlg.h.

Constructor & Destructor Documentation

cluster::BlurredClusteringAlg::BlurredClusteringAlg

(

fhicl::ParameterSet const &

pset

)

Definition at line 37 of file BlurredClusteringAlg.cxx.

  : fDebug{pset.get<bool>("Debug", false)}
  , fDetector{pset.get<std::string>("Detector", "dune35t")}
  , fBlurWire{pset.get<int>("BlurWire")}
  , fBlurTick{pset.get<int>("BlurTick")}
  , fSigmaWire{pset.get<double>("SigmaWire")}
  , fSigmaTick{pset.get<double>("SigmaTick")}
  , fMaxTickWidthBlur{pset.get<int>("MaxTickWidthBlur")}
  , fClusterWireDistance{pset.get<int>("ClusterWireDistance")}
  , fClusterTickDistance{pset.get<int>("ClusterTickDistance")}
  , fNeighboursThreshold{pset.get<unsigned int>("NeighboursThreshold")}
  , fMinNeighbours{pset.get<unsigned int>("MinNeighbours")}
  , fMinSize{pset.get<unsigned int>("MinSize")}
  , fMinSeed{pset.get<double>("MinSeed")}
  , fTimeThreshold{pset.get<double>("TimeThreshold")}
  , fChargeThreshold{pset.get<double>("ChargeThreshold")}
  , fKernelWidth{2 * fBlurWire + 1}
  , fKernelHeight{2 * fBlurTick * fMaxTickWidthBlur + 1}
  , fAllKernels{MakeKernels()}
{}

cluster::BlurredClusteringAlg::~BlurredClusteringAlg

(

)

Definition at line 58 of file BlurredClusteringAlg.cxx.

{
  if (fDebugCanvas) {
    std::string closeName = fDebugPDFName;
    closeName.append("]");
    fDebugCanvas->Print(closeName.c_str());
    delete fDebugCanvas;
  }
}

Member Function Documentation

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

Takes a vector of clusters (itself a vector of hits) and turns them into clusters using the initial hit selection.

Definition at line 110 of file BlurredClusteringAlg.cxx.

{
  // Loop through the clusters (each a vector of bins)
  for (auto const& bins : allClusterBins) {
    // Convert the clusters (vectors of bins) to hits in a vector of recob::Hits
    art::PtrVector<recob::Hit> clusHits = ConvertBinsToRecobHits(image, bins);

    mf::LogInfo("BlurredClustering") << "Cluster made from " << bins.size() << " bins, of which "
                                     << clusHits.size() << " were real hits";

    // Make sure the clusters are above the minimum cluster size
    if (clusHits.size() < fMinSize) {
      mf::LogVerbatim("BlurredClustering")
        << "Cluster of size " << clusHits.size()
        << " not saved since it is smaller than the minimum cluster size, set to " << fMinSize;
      continue;
    }

    clusters.push_back(clusHits);
  }
}

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

private

Converts a vector of bins into a hit selection - not all the hits in the bins vector are real hits.

Definition at line 616 of file BlurredClusteringAlg.cxx.

{
  // Create the vector of hits to output
  art::PtrVector<recob::Hit> hits;

  // Look through the hits in the cluster
  for (auto const bin : bins) {
    // Take each hit and convert it to a recob::Hit
    art::Ptr<recob::Hit> const hit = ConvertBinToRecobHit(image, bin);

    // If this hit was a real hit put it in the hit selection
    if (!hit.isNull()) hits.push_back(hit);
  }

  // Return the vector of hits to make cluster
  return hits;
}

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

private

Returns the charge stored in the global bin value.

Definition at line 653 of file BlurredClusteringAlg.cxx.

{
  int const x = bin % image.size();
  int const y = bin / image.size();
  return image.at(x).at(y);
}

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

private

Converts a bin into a recob::Hit (not all of these bins correspond to recob::Hits - some are fake hits created by the blurring)

Definition at line 636 of file BlurredClusteringAlg.cxx.

{
  int const wire = bin % image.size();
  int const tick = bin / image.size();
  return fHitMap[wire][tick];
}

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

Takes hit map and returns a 2D vector representing wire and tick, filled with the charge.

Definition at line 136 of file BlurredClusteringAlg.cxx.

{
  // Define the size of this particular plane -- dynamically to avoid huge histograms
  int lowerTick = readoutWindowSize, upperTick{}, lowerWire = fGeom->MaxWires(), upperWire{};
  using lar::to_element;
  using ranges::views::transform;
  for (auto const& hit : hits | transform(to_element)) {
    int histWire = GlobalWire(hit.WireID());
    if (hit.PeakTime() < lowerTick) lowerTick = hit.PeakTime();
    if (hit.PeakTime() > upperTick) upperTick = hit.PeakTime();
    if (histWire < lowerWire) lowerWire = histWire;
    if (histWire > upperWire) upperWire = histWire;
  }
  fLowerTick = lowerTick - 20;
  fUpperTick = upperTick + 20;
  fLowerWire = lowerWire - 20;
  fUpperWire = upperWire + 20;

  // Use a map to keep a track of the real hits and their wire/ticks
  fHitMap.clear();
  fHitMap.resize(fUpperWire - fLowerWire,
                 std::vector<art::Ptr<recob::Hit>>(fUpperTick - fLowerTick));

  // Create a 2D vector
  std::vector<std::vector<double>> image(fUpperWire - fLowerWire,
                                         std::vector<double>(fUpperTick - fLowerTick));

  // Look through the hits
  for (auto const& hit : hits) {
    int const wire = GlobalWire(hit->WireID());
    auto const tick = static_cast<int>(hit->PeakTime());
    float const charge = hit->Integral();

    // Fill hit map and keep a note of all real hits for later
    if (charge > image.at(wire - fLowerWire).at(tick - fLowerTick)) {
      image.at(wire - fLowerWire).at(tick - fLowerTick) = charge;
      fHitMap[wire - fLowerWire][tick - fLowerTick] = hit;
    }
  }

  // Keep a note of dead wires
  fDeadWires = std::vector<bool>(fUpperWire - fLowerWire, false);
  geo::PlaneID const planeID = hits.front()->WireID().planeID();

  for (int wire = fLowerWire; wire < fUpperWire; ++wire) {
    raw::ChannelID_t const channel =
      fGeom->PlaneWireToChannel(planeID.Plane, wire, planeID.TPC, planeID.Cryostat);
    fDeadWires[wire - fLowerWire] = !fChanStatus.IsGood(channel);
  }

  return image;
}

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

private

Converts an xbin and a ybin to a global bin number.

Definition at line 645 of file BlurredClusteringAlg.cxx.

{
  return ybin * image.size() + xbin;
}

void cluster::BlurredClusteringAlg::CreateDebugPDF	(	int	run,
		int	subrun,
		int	event
	)

Create the PDF to save debug images.

Definition at line 69 of file BlurredClusteringAlg.cxx.

{
  if (!fDebugCanvas) {

    // Create the grayscale palette for the Z axis
    Double_t Red[2] = {1.00, 0.00};
    Double_t Green[2] = {1.00, 0.00};
    Double_t Blue[2] = {1.00, 0.00};
    Double_t Length[2] = {0.00, 1.00};
    TColor::CreateGradientColorTable(2, Length, Red, Green, Blue, 1000);
    gStyle->SetOptStat(110000);

    // Decide what to call this PDF
    std::ostringstream oss;
    oss << "BlurredImages_Run" << run << "_Subrun" << subrun;
    fDebugPDFName = oss.str();
    fDebugCanvas = new TCanvas(fDebugPDFName.c_str(), "Image canvas", 1000, 500);
    fDebugPDFName.append(".pdf");

    std::string openName = fDebugPDFName;
    openName.append("[");
    fDebugCanvas->Print(openName.c_str());
    fDebugCanvas->Divide(2, 2);
    fDebugCanvas->SetGrid();
  }

  // Clear the pads on the canvas
  for (int i = 1; i <= 4; ++i) {
    fDebugCanvas->GetPad(i)->Clear();
  }

  std::ostringstream oss;
  oss << "Event " << event;
  fDebugCanvas->cd(1);
  TLatex l;
  l.SetTextSize(0.15);
  l.DrawLatex(0.1, 0.1, oss.str().c_str());
  fDebugCanvas->Print(fDebugPDFName.c_str());
}

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

private

Count how many dead wires there are in the blurring region for a particular hit Returns a pair of counters representing how many dead wires there are below and above the hit respectively

Definition at line 662 of file BlurredClusteringAlg.cxx.

{
  auto deadWires = std::make_pair(0, 0);

  int const lower_bin = width / 2;
  int const upper_bin = (width + 1) / 2;

  auto const offset = wire_bin + fLowerWire;
  for (int wire = std::max(offset - lower_bin, fLowerWire);
       wire < std::min(offset + upper_bin, fUpperWire);
       ++wire) {
    if (!fDeadWires[wire - fLowerWire]) continue;

    if (wire < offset)
      ++deadWires.first;
    else if (wire > offset)
      ++deadWires.second;
  }

  return deadWires;
}

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

private

Dynamically find the blurring radii and Gaussian sigma in each dimension.

Definition at line 685 of file BlurredClusteringAlg.cxx.

{
  // Calculate least squares slope
  double nhits{}, sumx{}, sumy{}, sumx2{}, sumxy{};
  for (unsigned int wireIt = 0; wireIt < fHitMap.size(); ++wireIt) {
    for (unsigned int tickIt = 0; tickIt < fHitMap.at(wireIt).size(); ++tickIt) {
      if (fHitMap[wireIt][tickIt].isNull()) continue;
      ++nhits;
      int const x = wireIt + fLowerWire;
      int const y = tickIt + fLowerTick;
      sumx += x;
      sumy += y;
      sumx2 += x * x;
      sumxy += x * y;
    }
  }
  double const gradient = (nhits * sumxy - sumx * sumy) / (nhits * sumx2 - sumx * sumx);

  // Get the rough unit vector for the trajectories, making sure to
  // catch the vertical gradient.
  auto const unit = std::isnan(gradient) ? TVector2{0, 1} : TVector2{1, gradient}.Unit();

  // Use this direction to scale the blurring radii and Gaussian sigma
  int const blur_wire = std::max(std::abs(std::round(fBlurWire * unit.X())), 1.);
  int const blur_tick = std::max(std::abs(std::round(fBlurTick * unit.Y())), 1.);

  int const sigma_wire = std::max(std::abs(std::round(fSigmaWire * unit.X())), 1.);
  int const sigma_tick = std::max(std::abs(std::round(fSigmaTick * unit.Y())), 1.);
  return {{blur_wire, blur_tick, sigma_wire, sigma_tick}};
}

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

Find clusters in the histogram.

Definition at line 192 of file BlurredClusteringAlg.cxx.

{
  // Size of image in x and y
  int const nbinsx = blurred.size();
  int const nbinsy = blurred.at(0).size();
  int const nbins = nbinsx * nbinsy;

  // Vectors to hold hit information
  std::vector<bool> used(nbins);
  std::vector<std::pair<double, int>> values;

  // Place the bin number and contents as a pair in the values vector
  for (int xbin = 0; xbin < nbinsx; ++xbin) {
    for (int ybin = 0; ybin < nbinsy; ++ybin) {
      int const bin = ConvertWireTickToBin(blurred, xbin, ybin);
      values.emplace_back(ConvertBinToCharge(blurred, bin), bin);
    }
  }

  // Sort the values into charge order
  std::sort(values.rbegin(), values.rend());

  // Count the number of iterations of the cluster forming loop (== number of clusters)
  int niter = 0;

  // Clustering loops
  // First loop - considers highest charge hits in decreasing order, and puts them in a new cluster if they aren't already clustered (makes new cluster every iteration)
  // Second loop - looks at the direct neighbours of this seed and clusters to this if above charge/time thresholds. Runs recursively over all hits in cluster (inc. new ones)
  while (true) {

    // Start a new cluster each time loop is executed
    std::vector<int> cluster;
    std::vector<double> times;

    // Get the highest charge bin (go no further if below seed threshold)
    if (double const blurred_binval = values[niter].first; blurred_binval < fMinSeed) break;

    // Iterate through the bins from highest charge down
    int const bin = values[niter++].second;

    // Put this bin in used if not already there
    if (used[bin]) continue;
    used[bin] = true;

    // Start a new cluster
    cluster.push_back(bin);

    // Get the time of this hit
    if (double const time = GetTimeOfBin(blurred, bin); time > 0) times.push_back(time);

    // Now cluster neighbouring hits to this seed
    while (true) {

      bool added_cluster{false};

      for (unsigned int clusBin = 0; clusBin < cluster.size(); ++clusBin) {

        // Get x and y values for bin (c++ returns a%b = a if a<b)
        int const binx = cluster[clusBin] % nbinsx;
        int const biny = ((cluster[clusBin] - binx) / nbinsx) % nbinsy;

        // Look for hits in the neighbouring x/y bins
        for (int x = binx - fClusterWireDistance; x <= binx + fClusterWireDistance; x++) {
          if (x >= nbinsx or x < 0) continue;
          for (int y = biny - fClusterTickDistance; y <= biny + fClusterTickDistance; y++) {
            if (y >= nbinsy or y < 0) continue;
            if (x == binx and y == biny) continue;

            // Get this bin
            auto const bin = ConvertWireTickToBin(blurred, x, y);
            if (bin >= nbinsx * nbinsy or bin < 0) continue;
            if (used[bin]) continue;

            // Get the blurred value and time for this bin
            double const blurred_binval = ConvertBinToCharge(blurred, bin);
            double const time =
              GetTimeOfBin(blurred, bin); // NB for 'fake' hits, time is defaulted to -10000

            // Check real hits pass time cut (ignores fake hits)
            if (time > 0 && times.size() > 0 && !PassesTimeCut(times, time)) continue;

            // Add to cluster if bin value is above threshold
            if (blurred_binval > fChargeThreshold) {
              used[bin] = true;
              cluster.push_back(bin);
              added_cluster = true;
              if (time > 0) { times.push_back(time); }
            } // End of adding blurred bin to cluster
          }
        } // End of looking at directly neighbouring bins

      } // End of looping over bins already in this cluster

      if (!added_cluster) break;

    } // End of adding hits to this cluster

    // Check this cluster is above minimum size
    if (cluster.size() < fMinSize) {
      for (auto const bin : cluster) {
        assert(bin >= 0);
        used[bin] = false;
      }
      continue;
    }

    // Fill in holes in the cluster
    for (unsigned int clusBin = 0; clusBin < cluster.size(); clusBin++) {

      // Looks at directly neighbouring bins (and not itself)
      for (int x = -1; x <= 1; ++x) {
        for (int y = -1; y <= 1; ++y) {
          if (x == 0 && y == 0) continue;

          // Look at neighbouring bins to the clustered bin which are inside the cluster
          int neighbouringBin = cluster[clusBin] + x + (y * nbinsx);
          if (neighbouringBin < nbinsx || neighbouringBin % nbinsx == 0 ||
              neighbouringBin % nbinsx == nbinsx - 1 || neighbouringBin >= nbinsx * (nbinsy - 1))
            continue;

          double const time = GetTimeOfBin(blurred, neighbouringBin);

          // If not already clustered and passes neighbour/time thresholds, add to cluster
          if (!used[neighbouringBin] &&
              (NumNeighbours(nbinsx, used, neighbouringBin) > fNeighboursThreshold) &&
              PassesTimeCut(times, time)) {
            used[neighbouringBin] = true;
            cluster.push_back(neighbouringBin);

            if (time > 0) { times.push_back(time); }
          } // End of clustering neighbouring bin
        }
      } // End of looping over neighbouring bins

    } // End of looping over bins already in cluster

    mf::LogVerbatim("Blurred Clustering")
      << "Size of cluster after filling in holes: " << cluster.size();

    // Remove peninsulas - hits which have too few neighbouring hits in the cluster (defined by fMinNeighbours)
    while (true) {
      bool removed_cluster{false};

      // Loop over all the bins in the cluster
      for (int clusBin = cluster.size() - 1; clusBin >= 0; clusBin--) {
        auto const bin = cluster[clusBin];

        // If bin is in cluster ignore
        if (bin < nbinsx || bin % nbinsx == 0 || bin % nbinsx == nbinsx - 1 ||
            bin >= nbinsx * (nbinsy - 1))
          continue;

        // Remove hit if it has too few neighbouring hits
        if (NumNeighbours(nbinsx, used, bin) < fMinNeighbours) {
          used[bin] = false;
          removed_cluster = true;
          cluster.erase(cluster.begin() + clusBin);
        }
      }

      if (!removed_cluster) break;
    }

    mf::LogVerbatim("Blurred Clustering")
      << "Size of cluster after removing peninsulas: " << cluster.size();

    // Disregard cluster if not of minimum size
    if (cluster.size() < fMinSize) {
      for (auto const bin : cluster) {
        assert(bin >= 0);
        used[bin] = false;
      }
      continue;
    }

    // Put this cluster in the vector of clusters
    allcluster.push_back(cluster);

  } // End loop over this cluster

  // Return the number of clusters found in this hit map
  return allcluster.size();
}

std::vector< std::vector< double > > cluster::BlurredClusteringAlg::GaussianBlur

(

std::vector< std::vector< double >> const &

image

)

const

Applies Gaussian blur to image.

Definition at line 433 of file BlurredClusteringAlg.cxx.

{
  if (fSigmaWire == 0 and fSigmaTick == 0) return image;

  auto const [blur_wire, blur_tick, sigma_wire, sigma_tick] = FindBlurringParameters();

  // Convolve the Gaussian
  int width = 2 * blur_wire + 1;
  int height = 2 * blur_tick + 1;
  int nbinsx = image.size();
  int nbinsy = image.at(0).size();

  // Blurred histogram and normalisation for each bin
  std::vector<std::vector<double>> copy(nbinsx, std::vector<double>(nbinsy, 0));

  // Loop through all the bins in the histogram to blur
  for (int x = 0; x < nbinsx; ++x) {
    for (int y = 0; y < nbinsy; ++y) {

      if (image[x][y] == 0) continue;

      // Scale the tick blurring based on the width of the hit
      int tick_scale =
        std::sqrt(cet::square(fHitMap[x][y]->RMS()) + cet::square(sigma_tick)) / (double)sigma_tick;
      tick_scale = std::max(std::min(tick_scale, fMaxTickWidthBlur), 1);
      auto const& correct_kernel = fAllKernels[sigma_wire][sigma_tick * tick_scale];

      // Find any dead wires in the potential blurring region
      auto const [lower_bin_dead, upper_bin_dead] = DeadWireCount(x, width);

      // Note of how many dead wires we have passed whilst blurring in the wire direction
      // If blurring below the seed hit, need to keep a note of how many dead wires to come
      // If blurring above, need to keep a note of how many dead wires have passed
      auto dead_wires_passed{lower_bin_dead};

      // Loop over the blurring region around this hit
      for (int blurx = -(width / 2 + lower_bin_dead); blurx < (width + 1) / 2 + upper_bin_dead;
           ++blurx) {
        if (x + blurx < 0) continue;
        for (int blury = -height / 2 * tick_scale;
             blury < ((((height + 1) / 2) - 1) * tick_scale) + 1;
             ++blury) {
          if (blurx < 0 and fDeadWires[x + blurx]) dead_wires_passed -= 1;

          // Smear the charge of this hit
          double const weight = correct_kernel[fKernelWidth * (fKernelHeight / 2 + blury) +
                                               (fKernelWidth / 2 + (blurx - dead_wires_passed))];
          if (x + blurx >= 0 and x + blurx < nbinsx and y + blury >= 0 and y + blury < nbinsy)
            copy[x + blurx][y + blury] += weight * image[x][y];

          if (blurx > 0 and fDeadWires[x + blurx]) dead_wires_passed += 1;
        }
      } // blurring region
    }
  } // hits to blur

  // HAVE REMOVED NOMALISATION CODE
  // WHEN USING DIFFERENT KERNELS, THERE'S NO EASY WAY OF DOING THIS...
  // RECONSIDER...

  // Return the blurred histogram
  return copy;
}

unsigned int cluster::BlurredClusteringAlg::GetMinSize ( ) const

inlinenoexcept

Minimum size of cluster to save.

Definition at line 80 of file BlurredClusteringAlg.h.

  {
    return fMinSize;
  }

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

private

Returns the hit time of a hit in a particular bin.

Definition at line 717 of file BlurredClusteringAlg.cxx.

{
  auto const hit = ConvertBinToRecobHit(image, bin);
  return hit.isNull() ? -10000. : hit->PeakTime();
}

int cluster::BlurredClusteringAlg::GlobalWire

(

geo::WireID const &

wireID

)

const

Find the global wire position.

Definition at line 378 of file BlurredClusteringAlg.cxx.

{
  double globalWire = -999;

  // Induction
  if (fGeom->SignalType(wireID) == geo::kInduction) {
    double wireCentre[3];
    fGeom->WireIDToWireGeo(wireID).GetCenter(wireCentre);
    if (wireID.TPC % 2 == 0)
      globalWire =
        fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 0, wireID.Cryostat);
    else
      globalWire =
        fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 1, wireID.Cryostat);
  }

  // Collection
  else {
    // FOR COLLECTION WIRES, HARD CODE THE GEOMETRY FOR GIVEN DETECTORS
    // THIS _SHOULD_ BE TEMPORARY. GLOBAL WIRE SUPPORT IS BEING ADDED TO THE LARSOFT GEOMETRY AND SHOULD BE AVAILABLE SOON
    if (fDetector == "dune35t") {
      unsigned int nwires = fGeom->Nwires(wireID.Plane, 0, wireID.Cryostat);
      if (wireID.TPC == 0 or wireID.TPC == 1)
        globalWire = wireID.Wire;
      else if (wireID.TPC == 2 or wireID.TPC == 3 or wireID.TPC == 4 or wireID.TPC == 5)
        globalWire = nwires + wireID.Wire;
      else if (wireID.TPC == 6 or wireID.TPC == 7)
        globalWire = (2 * nwires) + wireID.Wire;
      else
        mf::LogError("BlurredClusterAlg")
          << "Error when trying to find a global induction plane coordinate for TPC " << wireID.TPC
          << " (geometry " << fDetector << ")";
    }
    else if (fDetector == "dune10kt") {
      unsigned int nwires = fGeom->Nwires(wireID.Plane, 0, wireID.Cryostat);
      // Detector geometry has four TPCs, two on top of each other, repeated along z...
      int block = wireID.TPC / 4;
      globalWire = (nwires * block) + wireID.Wire;
    }
    else {
      double wireCentre[3];
      fGeom->WireIDToWireGeo(wireID).GetCenter(wireCentre);
      if (wireID.TPC % 2 == 0)
        globalWire =
          fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 0, wireID.Cryostat);
      else
        globalWire =
          fGeom->WireCoordinate(wireCentre[1], wireCentre[2], wireID.Plane, 1, wireID.Cryostat);
    }
  }

  return std::round(globalWire);
}

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

Converts a 2D vector in a histogram for the debug pdf.

Definition at line 498 of file BlurredClusteringAlg.cxx.

{
  auto hist = new TH2F(name,
                       name,
                       fUpperWire - fLowerWire,
                       fLowerWire - 0.5,
                       fUpperWire - 0.5,
                       fUpperTick - fLowerTick,
                       fLowerTick - 0.5,
                       fUpperTick - 0.5);
  hist->SetXTitle("Wire number");
  hist->SetYTitle("Tick number");
  hist->SetZTitle("Charge");

  for (unsigned int imageWireIt = 0; imageWireIt < image.size(); ++imageWireIt) {
    int const wire = imageWireIt + fLowerWire;
    for (unsigned int imageTickIt = 0; imageTickIt < image.at(imageWireIt).size(); ++imageTickIt) {
      int const tick = imageTickIt + fLowerTick;
      hist->Fill(wire, tick, image.at(imageWireIt).at(imageTickIt));
    }
  }

  return hist;
}

std::vector< std::vector< std::vector< double > > > cluster::BlurredClusteringAlg::MakeKernels ( ) const

private

Makes all the kernels which could be required given the tuned parameters.

Definition at line 725 of file BlurredClusteringAlg.cxx.

{
  // Kernel size is the largest possible given the hit width rescaling
  std::vector<std::vector<std::vector<double>>> allKernels(
    fSigmaWire + 1,
    std::vector<std::vector<double>>(fSigmaTick * fMaxTickWidthBlur + 1,
                                     std::vector<double>(fKernelWidth * fKernelHeight)));

  // Ranges of kernels to make
  // Complete range of sigmas possible after dynamic fixing and hit width convolution
  for (int sigma_wire = 1; sigma_wire <= fSigmaWire; ++sigma_wire) {
    for (int sigma_tick = 1; sigma_tick <= fSigmaTick * fMaxTickWidthBlur; ++sigma_tick) {

      // New kernel
      std::vector<double> kernel(fKernelWidth * fKernelHeight, 0);

      // Smear out according to the blur radii in each direction
      for (int i = -fBlurWire; i <= fBlurWire; i++) {
        for (int j = -fBlurTick * fMaxTickWidthBlur; j <= fBlurTick * fMaxTickWidthBlur; j++) {

          // Fill kernel
          double const sig2i = 2. * sigma_wire * sigma_wire;
          double const sig2j = 2. * sigma_tick * sigma_tick;

          int const key = (fKernelWidth * (j + fBlurTick * fMaxTickWidthBlur)) + (i + fBlurWire);
          double const value = 1. / std::sqrt(sig2i * M_PI) * std::exp(-i * i / sig2i) * 1. /
                               std::sqrt(sig2j * M_PI) * std::exp(-j * j / sig2j);
          kernel.at(key) = value;
        }
      } // End loop over blurring region

      allKernels[sigma_wire][sigma_tick] = move(kernel);
    }
  }
  return allKernels;
}

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

private

Determines the number of clustered neighbours of a hit.

2D hists can be considered a string of bins - the equation to convert between them is [bin = x + (nbinsx * y)]

Definition at line 763 of file BlurredClusteringAlg.cxx.

{
  unsigned int neighbours = 0;

  // Loop over all directly neighbouring hits (not itself)
  for (int x = -1; x <= 1; x++) {
    for (int y = -1; y <= 1; y++) {
      if (x == 0 && y == 0) continue;

      // Determine bin
      int neighbouringBin =
        bin + x +
        (y *
         nbinsx); /// 2D hists can be considered a string of bins - the equation to convert between them is [bin = x + (nbinsx * y)]

      // If this bin is in the cluster, increase the neighbouring bin counter
      if (used.at(neighbouringBin)) neighbours++;
    }
  }

  // Return the number of neighbours in the cluster of a particular hit
  return neighbours;
}

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

private

Determine if a hit is within a time threshold of any other hits in a cluster.

Definition at line 790 of file BlurredClusteringAlg.cxx.

{
  for (auto const t : times) {
    if (std::abs(time - t) < fTimeThreshold) return true;
  }
  return false;
}

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

Save the images for debugging This version takes the final clusters and overlays on the hit map

Definition at line 525 of file BlurredClusteringAlg.cxx.

{
  // Make a vector of clusters
  std::vector<std::vector<int>> allClusterBins;

  for (auto const& cluster : allClusters) {
    if (cluster.empty()) continue;

    std::vector<int> clusterBins;

    for (auto const& hit : cluster) {
      unsigned int const wire = GlobalWire(hit->WireID());
      float const tick = hit->PeakTime();
      int bin = image->GetBin((wire - fLowerWire) + 1, (tick - fLowerTick) + 1);
      if (cluster.size() < fMinSize) bin *= -1;

      clusterBins.push_back(bin);
    }

    allClusterBins.push_back(clusterBins);
  }

  SaveImage(image, allClusterBins, pad, tpc, plane);
}

void cluster::BlurredClusteringAlg::SaveImage	(	TH2F *	image,
		int	pad,
		int	tpc,
		int	plane
	)

Save the images for debugging.

Definition at line 555 of file BlurredClusteringAlg.cxx.

{
  std::vector<std::vector<int>> allClusterBins;
  SaveImage(image, allClusterBins, pad, tpc, plane);
}

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

Save the images for debugging This version takes a vector of bins and overlays the relevant bins on the hit map

Definition at line 562 of file BlurredClusteringAlg.cxx.

{
  fDebugCanvas->cd(pad);
  std::string stage;

  switch (pad) {
  case 1: stage = "Stage 1: Unblurred"; break;
  case 2: stage = "Stage 2: Blurred"; break;
  case 3: stage = "Stage 3: Blurred with clusters overlaid"; break;
  case 4: stage = "Stage 4: Output clusters"; break;
  default: stage = "Unknown stage"; break;
  }

  std::stringstream title;
  title << stage << " -- TPC " << tpc << ", Plane " << plane; // << " (Event " << fEvent << ")";

  image->SetName(title.str().c_str());
  image->SetTitle(title.str().c_str());
  image->DrawCopy("colz");

  // Draw the clustered hits on the histograms
  int clusterNum = 2;
  for (auto const& bins : allClusterBins) {
    TMarker mark(0, 0, 20);
    mark.SetMarkerColor(clusterNum);
    mark.SetMarkerSize(0.1);

    for (auto bin : bins) {
      // Hit from a cluster that we aren't going to save
      if (bin < 0) {
        bin *= -1;
        mark.SetMarkerStyle(24);
      }

      int wire, tick, z;
      image->GetBinXYZ(bin, wire, tick, z);
      mark.DrawMarker(wire + fLowerWire - 1, tick + fLowerTick - 1);
      mark.SetMarkerStyle(20);
    }
  }

  if (pad == 4) {
    fDebugCanvas->Print(fDebugPDFName.c_str());
    fDebugCanvas->Clear("D");
  }
}

Member Data Documentation

std::vector<std::vector<std::vector<double> > > cluster::BlurredClusteringAlg::fAllKernels

private

Definition at line 161 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fBlurTick

private

Definition at line 146 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fBlurWire

private

Definition at line 145 of file BlurredClusteringAlg.h.

lariov::ChannelStatusProvider const& cluster::BlurredClusteringAlg::fChanStatus

private

Initial value:

{
    art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider()}

Definition at line 176 of file BlurredClusteringAlg.h.

double cluster::BlurredClusteringAlg::fChargeThreshold

private

Definition at line 157 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fClusterTickDistance

private

Definition at line 151 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fClusterWireDistance

private

Definition at line 150 of file BlurredClusteringAlg.h.

std::vector<bool> cluster::BlurredClusteringAlg::fDeadWires

private

Definition at line 165 of file BlurredClusteringAlg.h.

bool cluster::BlurredClusteringAlg::fDebug

private

Definition at line 141 of file BlurredClusteringAlg.h.

TCanvas* cluster::BlurredClusteringAlg::fDebugCanvas {nullptr}

private

Definition at line 171 of file BlurredClusteringAlg.h.

std::string cluster::BlurredClusteringAlg::fDebugPDFName {}

private

Definition at line 172 of file BlurredClusteringAlg.h.

std::string cluster::BlurredClusteringAlg::fDetector

private

Definition at line 142 of file BlurredClusteringAlg.h.

art::ServiceHandle<geo::Geometry const> cluster::BlurredClusteringAlg::fGeom

private

Definition at line 175 of file BlurredClusteringAlg.h.

std::vector<std::vector<art::Ptr<recob::Hit> > > cluster::BlurredClusteringAlg::fHitMap

private

Definition at line 164 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fKernelHeight

private

Definition at line 160 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fKernelWidth

private

Definition at line 160 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fLowerTick

private

Definition at line 167 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fLowerWire

private

Definition at line 168 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fMaxTickWidthBlur

private

Definition at line 149 of file BlurredClusteringAlg.h.

unsigned int cluster::BlurredClusteringAlg::fMinNeighbours

private

Definition at line 153 of file BlurredClusteringAlg.h.

double cluster::BlurredClusteringAlg::fMinSeed

private

Definition at line 155 of file BlurredClusteringAlg.h.

unsigned int cluster::BlurredClusteringAlg::fMinSize

private

Definition at line 154 of file BlurredClusteringAlg.h.

unsigned int cluster::BlurredClusteringAlg::fNeighboursThreshold

private

Definition at line 152 of file BlurredClusteringAlg.h.

double cluster::BlurredClusteringAlg::fSigmaTick

private

Definition at line 148 of file BlurredClusteringAlg.h.

double cluster::BlurredClusteringAlg::fSigmaWire

private

Definition at line 147 of file BlurredClusteringAlg.h.

double cluster::BlurredClusteringAlg::fTimeThreshold

private

Definition at line 156 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fUpperTick

private

Definition at line 167 of file BlurredClusteringAlg.h.

int cluster::BlurredClusteringAlg::fUpperWire

private

Definition at line 168 of file BlurredClusteringAlg.h.

---

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

• larreco/larreco/RecoAlg/BlurredClusteringAlg.h
• larreco/larreco/RecoAlg/BlurredClusteringAlg.cxx