LArSoft  v10_04_05
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
43 
44 // ROOT & C++ includes
45 #include "TH2F.h"
46 #include <map>
47 #include <string>
48 
49 namespace cluster {
50  class BlurredClustering;
51 }
52 
54 public:
55  explicit BlurredClustering(fhicl::ParameterSet const& pset);
56 
57 private:
58  void produce(art::Event& evt) override;
59 
62 
63  // Create instances of algorithm classes to perform the clustering
67 };
68 
70  : EDProducer{pset}
71  , fHitsModuleLabel{pset.get<std::string>("HitsModuleLabel")}
72  , fTrackModuleLabel{pset.get<std::string>("TrackModuleLabel")}
73  , fVertexModuleLabel{pset.get<std::string>("VertexModuleLabel")}
74  , fPFParticleModuleLabel{pset.get<std::string>("PFParticleModuleLabel")}
75  , fCreateDebugPDF{pset.get<bool>("CreateDebugPDF")}
76  , fMergeClusters{pset.get<bool>("MergeClusters")}
77  , fGlobalTPCRecon{pset.get<bool>("GlobalTPCRecon")}
78  , fShowerReconOnly{pset.get<bool>("ShowerReconOnly")}
79  , fBlurredClusteringAlg{pset.get<fhicl::ParameterSet>("BlurredClusterAlg")}
80  , fMergeClusterAlg{pset.get<fhicl::ParameterSet>("MergeClusterAlg")}
81  , fTrackShowerSeparationAlg{pset.get<fhicl::ParameterSet>("TrackShowerSeparationAlg")}
82 {
83  produces<std::vector<recob::Cluster>>();
84  produces<art::Assns<recob::Cluster, recob::Hit>>();
85 }
86 
88 {
89  // Create debug pdf to illustrate the blurring process
91 
92  // Output containers -- collection of clusters and associations
93  auto clusters = std::make_unique<std::vector<recob::Cluster>>();
94  auto associations = std::make_unique<art::Assns<recob::Cluster, recob::Hit>>();
95 
96  // Compute the cluster characteristics
97  // Just use default for now, but configuration will go here
99 
100  // Services
102  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
103  auto const detProp =
105  util::GeometryUtilities const gser{
106  *geom, art::ServiceHandle<geo::WireReadout const>()->Get(), clockData, detProp};
107  int const readoutWindowSize = detProp.ReadOutWindowSize();
108 
109  // Get the hits from the event
111  std::vector<art::Ptr<recob::Hit>> hits;
112  std::vector<art::Ptr<recob::Hit>> hitsToCluster;
113  if (evt.getByLabel(fHitsModuleLabel, hitCollection)) art::fill_ptr_vector(hits, hitCollection);
114 
115  if (fShowerReconOnly) {
116 
117  // Get the tracks from the event
118  art::Handle<std::vector<recob::Track>> trackCollection;
119  std::vector<art::Ptr<recob::Track>> tracks;
120  if (evt.getByLabel(fTrackModuleLabel, trackCollection))
121  art::fill_ptr_vector(tracks, trackCollection);
122 
123  // Get the space points from the event
124  art::Handle<std::vector<recob::SpacePoint>> spacePointCollection;
125  std::vector<art::Ptr<recob::SpacePoint>> spacePoints;
126  if (evt.getByLabel(fTrackModuleLabel, spacePointCollection))
127  art::fill_ptr_vector(spacePoints, spacePointCollection);
128 
129  // Get vertices from the event
130  art::Handle<std::vector<recob::Vertex>> vertexCollection;
131  std::vector<art::Ptr<recob::Vertex>> vertices;
132  if (evt.getByLabel(fVertexModuleLabel, vertexCollection))
133  art::fill_ptr_vector(vertices, vertexCollection);
134 
135  // Get pandora pfparticles and clusters from the event
136  art::Handle<std::vector<recob::PFParticle>> pfParticleCollection;
137  std::vector<art::Ptr<recob::PFParticle>> pfParticles;
138  if (evt.getByLabel(fPFParticleModuleLabel, pfParticleCollection))
139  art::fill_ptr_vector(pfParticles, pfParticleCollection);
140  art::Handle<std::vector<recob::Cluster>> clusterCollection;
141  evt.getByLabel(fPFParticleModuleLabel, clusterCollection);
142 
143  if (trackCollection.isValid()) {
144  art::FindManyP<recob::Hit> fmht(trackCollection, evt, fTrackModuleLabel);
145  art::FindManyP<recob::Track> fmth(hitCollection, evt, fTrackModuleLabel);
146  art::FindManyP<recob::SpacePoint> fmspt(trackCollection, evt, fTrackModuleLabel);
147  art::FindManyP<recob::Track> fmtsp(spacePointCollection, evt, fTrackModuleLabel);
149  evt.event(), hits, tracks, spacePoints, fmht, fmth, fmspt, fmtsp);
150  }
151  }
152  else
153  hitsToCluster = hits;
154 
155  // Make a map between the planes and the hits on each
156  std::map<std::pair<int, int>, std::vector<art::Ptr<recob::Hit>>> planeToHits;
157  for (auto const& hitToCluster : hitsToCluster) {
158  auto const& wireID = hitToCluster->WireID();
159  auto const planeNo = wireID.Plane;
160  auto const tpc = fGlobalTPCRecon ? wireID.TPC % 2 : wireID.TPC;
161  planeToHits[std::make_pair(planeNo, tpc)].push_back(hitToCluster);
162  }
163 
164  // Loop over views
165  for (auto const& [plane, hits] : planeToHits) {
166  std::vector<art::PtrVector<recob::Hit>> finalClusters;
167 
168  // Implement the algorithm
169  if (hits.size() >= fBlurredClusteringAlg.GetMinSize()) {
170 
171  // Convert hit map to TH2 histogram and blur it
172  auto const image = fBlurredClusteringAlg.ConvertRecobHitsToVector(hits, readoutWindowSize);
173  auto const blurred = fBlurredClusteringAlg.GaussianBlur(image);
174 
175  // Find clusters in histogram
176  std::vector<std::vector<int>>
177  allClusterBins; // Vector of clusters (clusters are vectors of hits)
178  int numClusters = fBlurredClusteringAlg.FindClusters(blurred, allClusterBins);
179  mf::LogVerbatim("Blurred Clustering") << "Found " << numClusters << " clusters" << std::endl;
180 
181  // Create output clusters from the vector of clusters made in FindClusters
182  std::vector<art::PtrVector<recob::Hit>> planeClusters;
183  fBlurredClusteringAlg.ConvertBinsToClusters(image, allClusterBins, planeClusters);
184 
185  // Use the cluster merging algorithm
186  if (fMergeClusters) {
187  int numMergedClusters = fMergeClusterAlg.MergeClusters(planeClusters, finalClusters);
188  mf::LogVerbatim("Blurred Clustering")
189  << "After merging, there are " << numMergedClusters << " clusters" << std::endl;
190  }
191  else
192  finalClusters = planeClusters;
193 
194  // Make the debug PDF
195  if (fCreateDebugPDF) {
196  std::stringstream name;
197  name << "blurred_image";
198  TH2F* imageHist = fBlurredClusteringAlg.MakeHistogram(image, TString{name.str()});
199  name << "_convolved";
200  TH2F* blurredHist = fBlurredClusteringAlg.MakeHistogram(blurred, TString{name.str()});
201  auto const [planeNo, tpc] = plane;
202  fBlurredClusteringAlg.SaveImage(imageHist, 1, tpc, planeNo);
203  fBlurredClusteringAlg.SaveImage(blurredHist, 2, tpc, planeNo);
204  fBlurredClusteringAlg.SaveImage(blurredHist, allClusterBins, 3, tpc, planeNo);
205  fBlurredClusteringAlg.SaveImage(imageHist, finalClusters, 4, tpc, planeNo);
206  imageHist->Delete();
207  blurredHist->Delete();
208  }
209 
210  } // End min hits check
211 
212  // Make the output cluster objects
213  for (auto const& clusterHits : finalClusters) {
214  if (clusterHits.empty()) continue;
215 
216  // Get the start and end wires of the cluster
217  unsigned int const startWire =
218  fBlurredClusteringAlg.GlobalWire(clusterHits.front()->WireID());
219  unsigned int const endWire = fBlurredClusteringAlg.GlobalWire(clusterHits.back()->WireID());
220 
221  // Put cluster hits in the algorithm
222  ClusterParamAlgo.ImportHits(gser, clusterHits);
223 
224  // Create the recob::Cluster and place in the vector of clusters
225  ClusterCreator cluster(gser,
226  ClusterParamAlgo, // algo
227  float(startWire), // start_wire
228  0., // sigma_start_wire
229  clusterHits.front()->PeakTime(), // start_tick
230  clusterHits.front()->SigmaPeakTime(), // sigma_start_tick
231  float(endWire), // end_wire
232  0., // sigma_end_wire,
233  clusterHits.back()->PeakTime(), // end_tick
234  clusterHits.back()->SigmaPeakTime(), // sigma_end_tick
235  clusters->size(), // ID
236  clusterHits.front()->View(), // view
237  clusterHits.front()->WireID().planeID(), // plane
238  recob::Cluster::Sentry // sentry
239  );
240  clusters->emplace_back(cluster.move());
241 
242  // Associate the hits to this cluster
243  util::CreateAssn(evt, *clusters, clusterHits, *associations);
244  } // End loop over all clusters
245  }
246 
247  evt.put(std::move(clusters));
248  evt.put(std::move(associations));
249 }
250 
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
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
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.
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
Provides recob::Track data product.
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.