LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
SmallClusterFinderAlg.cxx
Go to the documentation of this file.
1 //
3 // \file SmallClusterFinderAlg.cxx
4 //
5 // \author corey.adams@yale.edu
6 //
7 // This algorithm is designed to find small clusters that could correspond to gammas
8 // or low energy electrons.
9 //
10 /* There are two parameters that matter from the fcl file:
11  fNHitsInClust is the number of hits that should be in these small clusters
12  ^-- Gamma events seem to rarely have more than 4 hits in the cluster
13  ^-- SN events are unclear. Should this even be used for SN?
14  fRadiusSizePar is the distance (in cm) between the small clusters and any other hits.
15 
16  This algorithm sorts the hits by plane, and then looks at each hit individually. If
17  there is a hit within RadiusSizePar, it gets added to a local list. All other hits
18  are ignored. Then, if the number of hits that got added to the local list is greater
19  then NHitsInClust, the original hit is ignored. If it's less, the original hit is
20  presumed to be part of a very small (or single hit) cluster. So its added to the list
21  of hits in the small cluster.
22 
23  All of the small clusters are then split apart into groups in the way you would expect.
24  Each cluster is assigned an ID number to distinguish it, and the hits that aren't
25  identified as small clusters all end up in the "leftover" cluster. The numbering scheme
26  is ID = 100*iPlane + Cluster on that plane, and the leftover hits are the first (0th)
27  cluster written out.
28 
29  -Corey
30 */
31 //
32 //
34 
35 #include <iostream>
36 
37 // Framework includes
40 #include "fhiclcpp/ParameterSet.h"
41 
42 // LArSoft includes
52 
53 // ***************** //
54 
55 //constructor with parameters:
57 {
58  fNPlanes = geom->Nplanes();
59  fRadiusSizePar = pset.get<double>("RadiusSizePar");
60  fNHitsInClust = pset.get<double>("NHitsInClust");
61  verbose = pset.get<bool>("Verbose");
63 }
64 
65 // ***************** //
66 // This method actually makes the clusters.
68  util::GeometryUtilities const& gser,
69  detinfo::DetectorClocksData const& clockData,
70  detinfo::DetectorPropertiesData const& detProp,
72 {
74  fDriftVelocity = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
76  fTimeTick = sampling_rate(clockData) / 1000.;
80 
82 
83  // Catch the case were there are no hits in the event:
84  if (allHits.size() == 0) {
85  if (verbose) std::cout << " no hits received! exiting " << std::endl;
86  return;
87  }
88 
89  art::Ptr<recob::Hit> theHit;
90 
91  //sort the hits into hits by plane
92  for (std::vector<art::Ptr<recob::Hit>>::iterator HitIter = allHits.begin();
93  HitIter != allHits.end();
94  HitIter++) {
95  theHit = *HitIter;
96  unsigned int p(0), w(0), t(0), cs(0); //c=channel, p=plane, w=wire, but what is t?
97  GetPlaneAndTPC(*HitIter, p, cs, t, w); //Find out what plane this hit is on.
98  //add this hit to the list specific to this plane
99  hitlistbyplane[p].push_back(theHit);
100  } // End loop on hits.
101 
102  // Want to check that the wires are OK for each event, and by each wire.
103  // This could be done in HitFinder, but I'm doing it here because I want to do
104  // some things specifically for Small Clusters... Flags to look for: Does
105  // this
106  // wire have more hits that the other wires that also have hits? Significantly
107  // more? Does the RMS of this hit lie well below the peak of the hit?
108  // Significantly below? More? If a bad wire is found, remove all of
109  // it's from the hitlists so they aren't clustered.
110 
111  //make the refined hit list for each plane.
112  for (unsigned int ip = 0; ip < fNPlanes; ip++) {
114 
115  //Check that the lists are populated correctly:
116  if (verbose)
117  std::cout << "At plane " << ip << ", found " << hitlistrefined[ip].size() << " hits with "
118  << hitlistleftover[ip].size() << " leftover" << std::endl;
119  //add the leftover hits to the correct object:
120  }
121 
122  //Now we have all of the gammas. This is good!
123  //Want to split all of the gammas up into individual clusters.
124 
125  //I was going to use lists instead of vectors to do this, because lists are more efficient
126  //at insertion and removal at arbitrary places. But, the number of gammas is so small that
127  //it just doesn't seem worth it
128  //
129  //Also, larsoft is so slow on its own that i think the difference is undetectable
130 
131  // The method to split the gammas into individual clusters that I want is:
132  // Use the existing method below to find all the hits within a certain distance from a given hit
133  // Take ALL of those hits and write them to a cluster.
134  // Then, remove those hits from the list of gammas.
135 
136  // Repeat, until the list of gammas is empty.
137 
138  // There is an issue of knowing which hits to remove from the total gamma hit list.
139  // To solve this, I'm going to overload SelectLocalHitList to take a reference to a vector
140  // of ints as an argument. When a hit is added to the local hit list,
141  // it's index will be added to a vector that is returned by reference.
142  // It's important to remove these hits in reverse order. This is because removing a hit
143  // changes the index of all of the hits after it
144 
145  //Now we need to take the lists of small clusters and sort it into the individual spots
146  //going to end up with smallClustList[plane][iClust][Hit]
147  //loop over planes of hits:
148 
149  for (unsigned int iplane = 0; iplane < fNPlanes; iplane++) {
150 
151  if (hitlistrefined[iplane].size() == 0) continue;
152 
153  //write the rest of the gammas one by one, in clusters:
154  int i = 1;
155 
156  std::vector<art::Ptr<recob::Hit>> splittingVector = hitlistrefined[iplane];
157 
158  while (splittingVector.size() != 0) {
159 
160  // std::cout << "\nThe hits remaining to be spilt are:" << std::endl;
161  //for (unsigned int j = 0; j < splittingVector.size();j++){
162  // std::cout << *splittingVector[j] << std::endl;
163  //}
164 
165  //find the first small cluster of gammas:
166  std::vector<int> index;
167  std::vector<art::Ptr<recob::Hit>> thiscluster;
168  thiscluster.clear();
169  index.clear();
170 
171  //Just use the first hit in the list of gammas:
172 
173  art::Ptr<recob::Hit> theHit = splittingVector.front(); //grab a hit from the list
174 
175  double time = theHit->PeakTime();
176  unsigned int plane(0), cstat(0), tpc(0), wire(0);
177  GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
178 
179  SelectLocalHitlist(gser, splittingVector, thiscluster, wire, time, fRadiusSizePar, index);
180 
181  if (verbose)
182  std::cout << "Done writing " << thiscluster.size() << " hits to cluster with ID "
183  << plane * 100 + i << std::endl;
184  //make sure to add these hits to the object that stores them:
185  smallClustList[plane].push_back(thiscluster);
186 
187  //Lastly, remove the gammas just clustered from the refinded hit list
188  while (index.size() != 0) {
189  splittingVector.erase(splittingVector.begin() + (index.back()));
190  index.pop_back();
191  }
192  i++;
193  }
194  }
195 }
196 
197 // ************************************* //
198 // Clear and resize - exactly what it sounds like
200 {
201  smallClustList.clear();
202  hitlistbyplane.clear();
203  hitlistrefined.clear();
204  hitlistleftover.clear();
205  smallClustList.resize(fNPlanes);
206  hitlistbyplane.resize(fNPlanes);
207  hitlistrefined.resize(fNPlanes);
208  hitlistleftover.resize(fNPlanes);
209 }
210 
211 /*
212  This is the method takes a list of hits ("hitlist") and compares each hit to a 2D
213  location that is specified by a wire ("wire_start") and time ("time_start"). If the
214  hit being examined is farther away than a specified distance ("radlimit", in cm) then
215  the hit is excluded. If the hit is within that distance, it's added.
216 */
218  util::GeometryUtilities const& gser,
220  std::vector<art::Ptr<recob::Hit>>& hitlistlocal,
221  double wire_start,
222  double time_start,
223  double radlimit) const
224 {
225  //loop over the hits in "hitlist", which should contain the hits we're selecting from
226  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hitIter = hitlist.begin();
227  hitIter != hitlist.end();
228  hitIter++) {
229  art::Ptr<recob::Hit> theHit = (*hitIter);
230  double time = theHit->PeakTime();
231  unsigned int plane, cstat, tpc, wire;
232  GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
233  //we now know which wire and what time this hit occurred on
234 
235  //calculate linear distance from start point and orthogonal distance from axis
236  double linear_dist = gser.Get2DDistance(wire, time, wire_start, time_start);
237 
238  if (linear_dist < radlimit) hitlistlocal.push_back(theHit);
239  }
240 }
241 
242 // This method is identical to the other method by the same name except that
243 // it keeps track of the location of the hits selected. That is, index is filled with
244 // the indices of the selected hits in the hitlist vector (the input vector)
246  util::GeometryUtilities const& gser,
248  std::vector<art::Ptr<recob::Hit>>& hitlistlocal,
249  double wire_start,
250  double time_start,
251  double radlimit,
252  std::vector<int>& index) const
253 {
254  // loop over the hits in "hitlist", which should contain the hits we're selecting from
255  int i = 0; //i keeps track of the index of the hit.
256  for (std::vector<art::Ptr<recob::Hit>>::const_iterator hitIter = hitlist.begin();
257  hitIter != hitlist.end();
258  hitIter++) {
259  art::Ptr<recob::Hit> theHit = (*hitIter);
260  double time = theHit->PeakTime();
261  unsigned int plane, cstat, tpc, wire;
262  GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
263  //we now know which wire and what time this hit occurred on
264 
265  //calculate linear distance from start point and orthogonal distance from axis
266  double linear_dist = gser.Get2DDistance(wire, time, wire_start, time_start);
267 
268  if (linear_dist < radlimit) {
269  hitlistlocal.push_back(theHit);
270  index.push_back(i);
271  }
272  i++;
273  }
274 
275  //need to make sure the index array is in order. Sort it!
276  std::sort(index.begin(), index.end());
277 }
278 
279 // This method sorts this hitlist (which should only be on a single plane)
280 std::vector<art::Ptr<recob::Hit>> cluster::SmallClusterFinderAlg::CreateHighHitlist(
281  util::GeometryUtilities const& gser,
282  std::vector<art::Ptr<recob::Hit>> const& hitlist,
283  std::vector<art::Ptr<recob::Hit>>& leftovers) const
284 {
285 
286  std::vector<art::Ptr<recob::Hit>>
287  hitlist_total; //This is the final result, a list of hits that are small clusters
288 
289  std::vector<art::Ptr<recob::Hit>> hitlistlocal;
290 
291  for (unsigned int ix = 0; ix < hitlist.size(); ix++) {
292 
293  art::Ptr<recob::Hit> const& theHit = hitlist[ix]; //grab a hit from the list
294 
295  double time = theHit->PeakTime();
296  unsigned int plane(0), cstat(0), tpc(0), wire(0);
297  //std::cout << "The hit is " << (*theHit) << std::endl;
298  GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
299  //use the wire and time of this hit as a seed.
300  // ^^^^^ This could probably be optimized?
301 
302  //get ALL of the hits from hitlist that are within the distance fRadiusSizePar of the seed hit.
303  SelectLocalHitlist(gser, hitlist, hitlistlocal, (double)wire, time, fRadiusSizePar);
304 
305  if (hitlistlocal.size() < fNHitsInClust) {
306  hitlist_total.push_back(theHit); //Add this hit if there are less than fNHitsInClust nearby.
307  if (verbose)
308  std::cout << " adding hit @ w,t " << wire << " " << time << " on plane " << plane
309  << std::endl;
310  }
311  else {
312  //Add this hit to the leftover pile
313  leftovers.push_back(theHit);
314  }
315 
316  hitlistlocal.clear(); //clear the local hit list, and look at the next hit.
317  }
318 
319  /*
320  This method could definitely be optimized. It creates a local hit list for each particle,
321  while if there is a local hit list that is sufficiently separated from all others it should be OK to
322  add them all at once. This is a task for future coders!
323  */
324 
325  return hitlist_total;
326 }
327 
328 // ******************************* //
330  unsigned int& p, //plane
331  unsigned int& /*cs*/, //cryostat
332  unsigned int& t, //time
333  unsigned int& w) const //wire
334 {
336  unsigned int channel = a->Channel();
337  geom->ChannelToWire(channel);
338  p = a->WireID().Plane;
339  t = a->PeakTime();
340  w = a->WireID().Wire;
341 
342  return 0;
343 }
344 
345 //Function to return the small clusters by plane
346 std::vector<std::vector<art::Ptr<recob::Hit>>>
348 {
349  if (iPlane < fNPlanes) return smallClustList[iPlane];
350  return {};
351 }
352 
353 std::vector<art::Ptr<recob::Hit>> cluster::SmallClusterFinderAlg::GetLeftoversByPlane(
354  unsigned int iPlane)
355 {
356  if (iPlane < fNPlanes) return hitlistleftover[iPlane];
357  return {};
358 }
std::vector< art::Ptr< recob::Hit > > CreateHighHitlist(util::GeometryUtilities const &gser, std::vector< art::Ptr< recob::Hit >> const &hitlist, std::vector< art::Ptr< recob::Hit >> &hitlistleftover) const
double Get2DDistance(double wire1, double time1, double wire2, double time2) const
Utilities related to art service access.
std::vector< WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
std::vector< std::vector< art::Ptr< recob::Hit > > > hitlistbyplane
std::vector< std::vector< std::vector< art::Ptr< recob::Hit > > > > smallClustList
Declaration of signal hit object.
double Temperature() const
In kelvin.
void FindSmallClusters(util::GeometryUtilities const &gser, detinfo::DetectorClocksData const &dataClocks, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> allHits)
int GetPlaneAndTPC(art::Ptr< recob::Hit > a, unsigned int &plane, unsigned int &cryostat, unsigned int &time, unsigned int &wire) const
void SelectLocalHitlist(util::GeometryUtilities const &gser, std::vector< art::Ptr< recob::Hit >> hitlist, std::vector< art::Ptr< recob::Hit >> &hitlistlocal, double wire_start, double time_start, double radlimit) const
std::vector< std::vector< art::Ptr< recob::Hit > > > hitlistleftover
double Efield(unsigned int planegap=0) const
kV/cm
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::vector< art::Ptr< recob::Hit > > GetLeftoversByPlane(unsigned int iPlane)
T get(std::string const &key) const
Definition: ParameterSet.h:314
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
std::vector< std::vector< art::Ptr< recob::Hit > > > hitlistrefined
Definition of data types for geometry description.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
Contains all timing reference information for the detector.
std::vector< std::vector< art::Ptr< recob::Hit > > > GetSmallClustersByPlane(unsigned int iPlane)
art::ServiceHandle< geo::Geometry const > geom
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
SmallClusterFinderAlg(fhicl::ParameterSet const &pset)
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
art framework interface to geometry description