#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 {}

geo::WireReadoutGeom const *	fWireReadoutGeom

lariov::ChannelStatusProvider const &	fChanStatus

Detailed Description

Definition at line 50 of file BlurredClusteringAlg.h.

Constructor & Destructor Documentation

cluster::BlurredClusteringAlg::BlurredClusteringAlg ( fhicl::ParameterSet const & pset )

Definition at line 37 of file BlurredClusteringAlg.cxx.

References fAllKernels, fBlurTick, fBlurWire, fChargeThreshold, fClusterTickDistance, fClusterWireDistance, fDetector, fKernelHeight, fKernelWidth, fMaxTickWidthBlur, fMinNeighbours, fMinSeed, fMinSize, fNeighboursThreshold, fSigmaTick, fSigmaWire, fTimeThreshold, fhicl::ParameterSet::get(), and MakeKernels().

   : 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.

References fDebugCanvas, and fDebugPDFName.

 {
   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 108 of file BlurredClusteringAlg.cxx.

References ConvertBinsToRecobHits(), fMinSize, and art::PtrVector< T >::size().

Referenced by cluster::BlurredClustering::produce().

 {
   // 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 604 of file BlurredClusteringAlg.cxx.

References bin, ConvertBinToRecobHit(), hits(), art::Ptr< T >::isNull(), and art::PtrVector< T >::push_back().

Referenced by ConvertBinsToClusters().

 {
   // 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 641 of file BlurredClusteringAlg.cxx.

References x, and y.

Referenced by FindClusters().

 {
   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 624 of file BlurredClusteringAlg.cxx.

References fHitMap.

Referenced by ConvertBinsToRecobHits(), and GetTimeOfBin().

 {
   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 133 of file BlurredClusteringAlg.cxx.

References fChanStatus, fDeadWires, fHitMap, fLowerTick, fLowerWire, fUpperTick, fUpperWire, fWireReadoutGeom, GlobalWire(), hits(), geo::WireReadoutGeom::MaxWires(), geo::WireReadoutGeom::PlaneWireToChannel(), lar::to_element, and lar::dump::vector().

Referenced by cluster::BlurredClustering::produce().

 {
   // Define the size of this particular plane -- dynamically to avoid huge histograms
   int lowerTick = readoutWindowSize, upperTick{}, lowerWire = fWireReadoutGeom->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();
 
   for (int wire = fLowerWire; wire < fUpperWire; ++wire) {
     raw::ChannelID_t const channel =
       fWireReadoutGeom->PlaneWireToChannel(geo::WireID(planeID, wire));
     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 633 of file BlurredClusteringAlg.cxx.

Referenced by FindClusters().

 {
   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 68 of file BlurredClusteringAlg.cxx.

References fDebugCanvas, fDebugPDFName, and tca::Length().

Referenced by cluster::BlurredClustering::produce().

 {
   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 650 of file BlurredClusteringAlg.cxx.

References fDeadWires, fLowerWire, and fUpperWire.

Referenced by GaussianBlur().

 {
   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 673 of file BlurredClusteringAlg.cxx.

References util::abs(), fBlurTick, fBlurWire, fHitMap, fLowerTick, fLowerWire, fSigmaTick, fSigmaWire, x, and y.

Referenced by GaussianBlur().

 {
   // 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 189 of file BlurredClusteringAlg.cxx.

References bin, ConvertBinToCharge(), ConvertWireTickToBin(), fChargeThreshold, fClusterTickDistance, fClusterWireDistance, fMinNeighbours, fMinSeed, fMinSize, fNeighboursThreshold, GetTimeOfBin(), nbinsx, NumNeighbours(), PassesTimeCut(), util::values(), x, and y.

Referenced by cluster::BlurredClustering::produce().

 {
   // 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 421 of file BlurredClusteringAlg.cxx.

References DeadWireCount(), fAllKernels, fDeadWires, fHitMap, FindBlurringParameters(), fKernelHeight, fKernelWidth, fMaxTickWidthBlur, fSigmaTick, fSigmaWire, nbinsx, weight, x, and y.

Referenced by cluster::BlurredClustering::produce().

 {
   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.

References bin, tca::DeadWireCount(), and lar::dump::vector().

Referenced by cluster::BlurredClustering::produce().

80 { return fMinSize; }

cluster::BlurredClusteringAlg::fMinSize

unsigned int fMinSize

Definition: BlurredClusteringAlg.h:151

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 704 of file BlurredClusteringAlg.cxx.

References ConvertBinToRecobHit().

Referenced by FindClusters().

 {
   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 374 of file BlurredClusteringAlg.cxx.

References geo::CryostatID::Cryostat, fDetector, fWireReadoutGeom, geo::WireGeo::GetCenter(), geo::kInduction, geo::WireReadoutGeom::Nwires(), geo::WireReadoutGeom::Plane(), geo::PlaneID::Plane, geo::WireReadoutGeom::SignalType(), geo::TPCID::TPC, geo::WireReadoutGeom::Wire(), and geo::WireID::Wire.

Referenced by ConvertRecobHitsToVector(), cluster::BlurredClustering::produce(), and SaveImage().

 {
   double globalWire = -999;
 
   // Induction
   if (fWireReadoutGeom->SignalType(wireID) == geo::kInduction) {
     auto const wireCenter = fWireReadoutGeom->Wire(wireID).GetCenter();
     globalWire =
       fWireReadoutGeom->Plane(geo::PlaneID{wireID.Cryostat, wireID.TPC % 2, wireID.Plane})
         .WireCoordinate(wireCenter);
   }
 
   // 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
     geo::PlaneID const planeid{wireID.Cryostat, 0, wireID.Plane};
     if (fDetector == "dune35t") {
       unsigned int nwires = fWireReadoutGeom->Nwires(planeid);
       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 = fWireReadoutGeom->Nwires(planeid);
       // 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 {
       auto const wireCenter = fWireReadoutGeom->Wire(wireID).GetCenter();
       globalWire =
         fWireReadoutGeom->Plane(geo::PlaneID{wireID.Cryostat, wireID.TPC % 2, wireID.Plane})
           .WireCoordinate(wireCenter);
     }
   }
 
   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 486 of file BlurredClusteringAlg.cxx.

References fLowerTick, fLowerWire, fUpperTick, fUpperWire, and hist.

Referenced by cluster::BlurredClustering::produce().

 {
   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 711 of file BlurredClusteringAlg.cxx.

References fBlurTick, fBlurWire, fKernelHeight, fKernelWidth, fMaxTickWidthBlur, fSigmaTick, fSigmaWire, value, and lar::dump::vector().

Referenced by BlurredClusteringAlg().

 {
   // 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 748 of file BlurredClusteringAlg.cxx.

References nbinsx, x, and y.

Referenced by FindClusters().

 {
   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); 
 
       // 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 774 of file BlurredClusteringAlg.cxx.

References util::abs(), and fTimeThreshold.

Referenced by FindClusters().

 {
   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 512 of file BlurredClusteringAlg.cxx.

References bin, fLowerTick, fLowerWire, fMinSize, and GlobalWire().

Referenced by cluster::BlurredClustering::produce(), and SaveImage().

 {
   // 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 542 of file BlurredClusteringAlg.cxx.

References SaveImage().

 {
   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 551 of file BlurredClusteringAlg.cxx.

References bin, fDebugCanvas, fDebugPDFName, fLowerTick, fLowerWire, and z.

 {
   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");
   }
 }