LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
BlurredClustering_module.cc
Go to the documentation of this file.
1 // Class: BlurredClustering
3 // Module Type: producer
4 // File: BlurredClustering_module.cc
5 // Author: Mike Wallbank (m.wallbank@sheffield.ac.uk), May 2015
6 //
7 // Reconstructs showers by blurring the hit map image to introduce fake
8 // hits before clustering to make fuller and more complete clusters.
9 //
10 // See DUNE-DocDB 54 (public) for details.
12 
13 // Framework includes:
21 #include "fhiclcpp/ParameterSet.h"
23 
24 // LArSoft includes
42 
43 // ROOT & C++ includes
44 #include "TH2F.h"
45 #include <map>
46 #include <string>
47 
48 namespace cluster {
49  class BlurredClustering;
50 }
51 
53 public:
54  explicit BlurredClustering(fhicl::ParameterSet const& pset);
55 
56 private:
57  void produce(art::Event& evt) override;
58 
61 
62  // Create instances of algorithm classes to perform the clustering
66 };
67 
69  : EDProducer{pset}
70  , fHitsModuleLabel{pset.get<std::string>("HitsModuleLabel")}
71  , fTrackModuleLabel{pset.get<std::string>("TrackModuleLabel")}
72  , fVertexModuleLabel{pset.get<std::string>("VertexModuleLabel")}
73  , fPFParticleModuleLabel{pset.get<std::string>("PFParticleModuleLabel")}
74  , fCreateDebugPDF{pset.get<bool>("CreateDebugPDF")}
75  , fMergeClusters{pset.get<bool>("MergeClusters")}
76  , fGlobalTPCRecon{pset.get<bool>("GlobalTPCRecon")}
77  , fShowerReconOnly{pset.get<bool>("ShowerReconOnly")}
78  , fBlurredClusteringAlg{pset.get<fhicl::ParameterSet>("BlurredClusterAlg")}
79  , fMergeClusterAlg{pset.get<fhicl::ParameterSet>("MergeClusterAlg")}
80  , fTrackShowerSeparationAlg{pset.get<fhicl::ParameterSet>("TrackShowerSeparationAlg")}
81 {
82  produces<std::vector<recob::Cluster>>();
83  produces<art::Assns<recob::Cluster, recob::Hit>>();
84 }
85 
87 {
88  // Create debug pdf to illustrate the blurring process
90 
91  // Output containers -- collection of clusters and associations
92  auto clusters = std::make_unique<std::vector<recob::Cluster>>();
93  auto associations = std::make_unique<art::Assns<recob::Cluster, recob::Hit>>();
94 
95  // Compute the cluster characteristics
96  // Just use default for now, but configuration will go here
98 
99  // Services
101  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
102  auto const detProp =
104  util::GeometryUtilities const gser{*geom, clockData, detProp};
105  int const readoutWindowSize = detProp.ReadOutWindowSize();
106 
107  // Get the hits from the event
109  std::vector<art::Ptr<recob::Hit>> hits;
110  std::vector<art::Ptr<recob::Hit>> hitsToCluster;
111  if (evt.getByLabel(fHitsModuleLabel, hitCollection)) art::fill_ptr_vector(hits, hitCollection);
112 
113  if (fShowerReconOnly) {
114 
115  // Get the tracks from the event
116  art::Handle<std::vector<recob::Track>> trackCollection;
117  std::vector<art::Ptr<recob::Track>> tracks;
118  if (evt.getByLabel(fTrackModuleLabel, trackCollection))
119  art::fill_ptr_vector(tracks, trackCollection);
120 
121  // Get the space points from the event
122  art::Handle<std::vector<recob::SpacePoint>> spacePointCollection;
123  std::vector<art::Ptr<recob::SpacePoint>> spacePoints;
124  if (evt.getByLabel(fTrackModuleLabel, spacePointCollection))
125  art::fill_ptr_vector(spacePoints, spacePointCollection);
126 
127  // Get vertices from the event
128  art::Handle<std::vector<recob::Vertex>> vertexCollection;
129  std::vector<art::Ptr<recob::Vertex>> vertices;
130  if (evt.getByLabel(fVertexModuleLabel, vertexCollection))
131  art::fill_ptr_vector(vertices, vertexCollection);
132 
133  // Get pandora pfparticles and clusters from the event
134  art::Handle<std::vector<recob::PFParticle>> pfParticleCollection;
135  std::vector<art::Ptr<recob::PFParticle>> pfParticles;
136  if (evt.getByLabel(fPFParticleModuleLabel, pfParticleCollection))
137  art::fill_ptr_vector(pfParticles, pfParticleCollection);
138  art::Handle<std::vector<recob::Cluster>> clusterCollection;
139  evt.getByLabel(fPFParticleModuleLabel, clusterCollection);
140 
141  if (trackCollection.isValid()) {
142  art::FindManyP<recob::Hit> fmht(trackCollection, evt, fTrackModuleLabel);
143  art::FindManyP<recob::Track> fmth(hitCollection, evt, fTrackModuleLabel);
144  art::FindManyP<recob::SpacePoint> fmspt(trackCollection, evt, fTrackModuleLabel);
145  art::FindManyP<recob::Track> fmtsp(spacePointCollection, evt, fTrackModuleLabel);
147  evt.event(), hits, tracks, spacePoints, fmht, fmth, fmspt, fmtsp);
148  }
149  }
150  else
151  hitsToCluster = hits;
152 
153  // Make a map between the planes and the hits on each
154  std::map<std::pair<int, int>, std::vector<art::Ptr<recob::Hit>>> planeToHits;
155  for (auto const& hitToCluster : hitsToCluster) {
156  auto const& wireID = hitToCluster->WireID();
157  auto const planeNo = wireID.Plane;
158  auto const tpc = fGlobalTPCRecon ? wireID.TPC % 2 : wireID.TPC;
159  planeToHits[std::make_pair(planeNo, tpc)].push_back(hitToCluster);
160  }
161 
162  // Loop over views
163  for (auto const& [plane, hits] : planeToHits) {
164  std::vector<art::PtrVector<recob::Hit>> finalClusters;
165 
166  // Implement the algorithm
167  if (hits.size() >= fBlurredClusteringAlg.GetMinSize()) {
168 
169  // Convert hit map to TH2 histogram and blur it
170  auto const image = fBlurredClusteringAlg.ConvertRecobHitsToVector(hits, readoutWindowSize);
171  auto const blurred = fBlurredClusteringAlg.GaussianBlur(image);
172 
173  // Find clusters in histogram
174  std::vector<std::vector<int>>
175  allClusterBins; // Vector of clusters (clusters are vectors of hits)
176  int numClusters = fBlurredClusteringAlg.FindClusters(blurred, allClusterBins);
177  mf::LogVerbatim("Blurred Clustering") << "Found " << numClusters << " clusters" << std::endl;
178 
179  // Create output clusters from the vector of clusters made in FindClusters
180  std::vector<art::PtrVector<recob::Hit>> planeClusters;
181  fBlurredClusteringAlg.ConvertBinsToClusters(image, allClusterBins, planeClusters);
182 
183  // Use the cluster merging algorithm
184  if (fMergeClusters) {
185  int numMergedClusters = fMergeClusterAlg.MergeClusters(planeClusters, finalClusters);
186  mf::LogVerbatim("Blurred Clustering")
187  << "After merging, there are " << numMergedClusters << " clusters" << std::endl;
188  }
189  else
190  finalClusters = planeClusters;
191 
192  // Make the debug PDF
193  if (fCreateDebugPDF) {
194  std::stringstream name;
195  name << "blurred_image";
196  TH2F* imageHist = fBlurredClusteringAlg.MakeHistogram(image, TString{name.str()});
197  name << "_convolved";
198  TH2F* blurredHist = fBlurredClusteringAlg.MakeHistogram(blurred, TString{name.str()});
199  auto const [planeNo, tpc] = plane;
200  fBlurredClusteringAlg.SaveImage(imageHist, 1, tpc, planeNo);
201  fBlurredClusteringAlg.SaveImage(blurredHist, 2, tpc, planeNo);
202  fBlurredClusteringAlg.SaveImage(blurredHist, allClusterBins, 3, tpc, planeNo);
203  fBlurredClusteringAlg.SaveImage(imageHist, finalClusters, 4, tpc, planeNo);
204  imageHist->Delete();
205  blurredHist->Delete();
206  }
207 
208  } // End min hits check
209 
210  // Make the output cluster objects
211  for (auto const& clusterHits : finalClusters) {
212  if (clusterHits.empty()) continue;
213 
214  // Get the start and end wires of the cluster
215  unsigned int const startWire =
216  fBlurredClusteringAlg.GlobalWire(clusterHits.front()->WireID());
217  unsigned int const endWire = fBlurredClusteringAlg.GlobalWire(clusterHits.back()->WireID());
218 
219  // Put cluster hits in the algorithm
220  ClusterParamAlgo.ImportHits(gser, clusterHits);
221 
222  // Create the recob::Cluster and place in the vector of clusters
223  ClusterCreator cluster(gser,
224  ClusterParamAlgo, // algo
225  float(startWire), // start_wire
226  0., // sigma_start_wire
227  clusterHits.front()->PeakTime(), // start_tick
228  clusterHits.front()->SigmaPeakTime(), // sigma_start_tick
229  float(endWire), // end_wire
230  0., // sigma_end_wire,
231  clusterHits.back()->PeakTime(), // end_tick
232  clusterHits.back()->SigmaPeakTime(), // sigma_end_tick
233  clusters->size(), // ID
234  clusterHits.front()->View(), // view
235  clusterHits.front()->WireID().planeID(), // plane
236  recob::Cluster::Sentry // sentry
237  );
238  clusters->emplace_back(cluster.move());
239 
240  // Associate the hits to this cluster
241  util::CreateAssn(evt, *clusters, clusterHits, *associations);
242  } // End loop over all clusters
243  }
244 
245  evt.put(std::move(clusters));
246  evt.put(std::move(associations));
247 }
248 
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
SubRunNumber_t subRun() const
Definition: Event.cc:35
int MergeClusters(std::vector< art::PtrVector< recob::Hit >> const &planeClusters, std::vector< art::PtrVector< recob::Hit >> &clusters) const
TH2F * MakeHistogram(std::vector< std::vector< double >> const &image, TString name) const
Converts a 2D vector in a histogram for the debug pdf.
Class managing the creation of a new recob::Cluster object.
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
cluster::MergeClusterAlg fMergeClusterAlg
Cluster finding and building.
int FindClusters(std::vector< std::vector< double >> const &image, std::vector< std::vector< int >> &allcluster) const
Find clusters in the histogram.
static const SentryArgument_t Sentry
An instance of the sentry object.
Definition: Cluster.h:174
bool isValid() const noexcept
Definition: Handle.h:203
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
void CreateDebugPDF(int run, int subrun, int event)
Create the PDF to save debug images.
void hits()
Definition: readHits.C:15
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
BlurredClustering(fhicl::ParameterSet const &pset)
Helper functions to create a cluster.
Wrapper for ClusterParamsAlgBase objects to accept diverse input.
int GlobalWire(geo::WireID const &wireID) const
Find the global wire position.
Provides recob::Track data product.
Wrapper for ClusterParamsAlgBase objects to accept arbitrary input.
std::vector< std::vector< double > > GaussianBlur(std::vector< std::vector< double >> const &image) const
Applies Gaussian blur to image.
EventNumber_t event() const
Definition: Event.cc:41
void SaveImage(TH2F *image, std::vector< art::PtrVector< recob::Hit >> const &allClusters, int pad, int tpc, int plane)
Declaration of cluster object.
void produce(art::Event &evt) override
shower::TrackShowerSeparationAlg fTrackShowerSeparationAlg
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
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 h...
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...
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
cluster::BlurredClusteringAlg fBlurredClusteringAlg
Interface to class computing cluster parameters.
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
std::vector< art::Ptr< recob::Hit > > SelectShowerHits(int event, const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< art::Ptr< recob::Track >> &tracks, const std::vector< art::Ptr< recob::SpacePoint >> &spacePoints, const art::FindManyP< recob::Hit > &fmht, const art::FindManyP< recob::Track > &fmth, const art::FindManyP< recob::SpacePoint > &fmspt, const art::FindManyP< recob::Track > &fmtsp) const
RunNumber_t run() const
Definition: Event.cc:29
art framework interface to geometry description
unsigned int GetMinSize() const noexcept
Minimum size of cluster to save.