LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
SmallClusterFinder_module.cc
Go to the documentation of this file.
1 //
3 // \file SmallClusterFinder.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  All of the heavy lifting is done is the SmallClusterFinderAlg.
30  This module remains responsible for interacting with the framework, getting hits,
31  writing clusters, etc.
32 
33  Thanks to Andrzej for the basic alg.
34 
35  -Corey
36 */
37 //
39 
40 #include <iostream>
41 
42 // include the proper bit of the framework
47 
48 //All the larsoft goodies:
62 
63 namespace cluster {
64 
66  public:
67  explicit SmallClusterFinder(fhicl::ParameterSet const& pset);
68 
69  private:
70  void beginJob() override;
71  void produce(art::Event& evt) override;
77 
78  //input parameters
79  std::string fHitFinderModuleLabel;
80  bool verbose;
81  // Remember, this is the *small* cluster finder
82 
83  SmallClusterFinderAlg fSmallClusterFinderAlg; //object that actually finds and sorts gammas
84 
85  unsigned int fNPlanes; // number of planes
86 
88  unsigned int& p,
89  unsigned int& cs,
90  unsigned int& t,
91  unsigned int& w);
92 
93  }; // class SmallAngleFinder
94 
95 }
96 namespace cluster {
98  : EDProducer{pset}, fSmallClusterFinderAlg(pset.get<fhicl::ParameterSet>("smallClustAlg"))
99  {
100  fHitFinderModuleLabel = pset.get<std::string>("HitFinderModuleLabel");
101  verbose = pset.get<bool>("Verbose");
102 
103  produces<std::vector<recob::Cluster>>(); //This code makes clusters
104  produces<art::Assns<recob::Cluster, recob::Hit>>(); //Matches clusters with hits
105  }
106 
107  // ***************** //
109  {
110  // this will not change on a run per run basis.
111  fNPlanes = wireReadoutGeom.Nplanes(); //get the number of planes in the ODTPC
112  }
113 
114  // ***************** //
115  // This method actually makes the clusters.
117  {
120  //Get the hits for this event:
122  evt.getByLabel(fHitFinderModuleLabel, HitListHandle);
123 
124  //A vector to hold hits, not yet filled:
125  std::vector<art::Ptr<recob::Hit>> hitlist;
126 
127  //How many hits in this event? Tell user:
128  if (verbose)
129  std::cout << " ++++ Hitsreceived received " << HitListHandle->size() << " +++++ "
130  << std::endl;
131  //Catch the case were there are no hits in the event:
132  if (HitListHandle->size() == 0) {
133  if (verbose) std::cout << "No hits received! Exiting." << std::endl;
134  return;
135  }
136  hitlist.resize(HitListHandle->size());
137 
138  //wrap the hits in art::Ptrs to pass to the Alg
139  for (unsigned int iHit = 0; iHit < hitlist.size(); iHit++) {
140  hitlist[iHit] = art::Ptr<recob::Hit>(HitListHandle, iHit);
141  }
142 
143  // Now run the alg to find the gammas:
144  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
145  auto const detProp =
147  util::GeometryUtilities const gser{*geom, wireReadoutGeom, clockData, detProp};
148  fSmallClusterFinderAlg.FindSmallClusters(gser, clockData, detProp, hitlist);
149 
150  // make an art::PtrVector of the clusters
151  auto SmallClusterFinder = std::make_unique<std::vector<recob::Cluster>>();
152  auto assn = std::make_unique<art::Assns<recob::Cluster, recob::Hit>>();
153 
154  // prepare the algorithm to compute the cluster characteristics;
155  // we use the "standard" one here; configuration would happen here,
156  // but we are using the default configuration for that algorithm
158 
159  for (unsigned int iplane = 0; iplane < fNPlanes; iplane++) {
160 
161  auto const leftoverHits = fSmallClusterFinderAlg.GetLeftoversByPlane(iplane);
162 
163  //write the leftover hits as a cluster:
164  if (leftoverHits.size() != 0) {
165  // pick some information from the first hit
166  geo::PlaneID planeID = leftoverHits.front()->WireID().planeID();
167  if (verbose)
168  std::cout << "Writing leftover hits to cluster ID: " << iplane * 100 << std::endl;
169 
170  ClusterParamAlgo.ImportHits(gser, leftoverHits);
171 
172  // create the recob::Cluster directly in the vector
173  ClusterCreator leftover(
174  gser,
175  ClusterParamAlgo, // algo
176  0., // start_wire
177  0., // sigma_start_wire
178  0., // start_tick
179  0., // sigma_start_tick
180  0., // end_wire
181  0., // sigma_end_wire,
182  0., // end_tick
183  0., // sigma_end_tick
184  iplane * 100, // ID
185  wireReadoutGeom.Plane(geo::PlaneID{planeID.asTPCID(), iplane}).View(),
186  planeID, // plane
187  recob::Cluster::Sentry // sentry
188  );
189 
190  SmallClusterFinder->emplace_back(leftover.move());
191 
192  util::CreateAssn(evt, *SmallClusterFinder, leftoverHits, *assn);
193  } //leftovers are written for this plane, if they exist.
194 
195  auto const smallClusters = fSmallClusterFinderAlg.GetSmallClustersByPlane(iplane);
196  for (unsigned int i = 0; i < smallClusters.size(); i++) {
197  // pick some information from the first hit
198  geo::PlaneID planeID; // invalid by default
199  if (!smallClusters.empty()) planeID = smallClusters[i].front()->WireID().planeID();
200 
201  ClusterParamAlgo.ImportHits(gser, smallClusters[i]);
202 
203  // create the recob::Cluster directly in the vector
204  ClusterCreator clust(gser,
205  ClusterParamAlgo, // algo
206  0., // start_wire
207  0., // sigma_start_wire
208  0., // start_tick
209  0., // sigma_start_tick
210  0., // end_wire
211  0., // sigma_end_wire,
212  0., // end_tick
213  0., // sigma_end_tick
214  iplane * 100 + i + 1, // ID
215  wireReadoutGeom.Plane(geo::PlaneID{planeID.asTPCID(), iplane}).View(),
216  planeID, // plane
217  recob::Cluster::Sentry // sentry
218  );
219 
220  SmallClusterFinder->emplace_back(clust.move());
221  // associate the hits to this cluster
222  util::CreateAssn(evt, *SmallClusterFinder, smallClusters[i], *assn);
223  }
224  }
225 
226  evt.put(std::move(SmallClusterFinder));
227  evt.put(std::move(assn));
228  } //end produce
229 
231 }
Class managing the creation of a new recob::Cluster object.
void produce(art::Event &evt) override
SmallClusterFinderAlg fSmallClusterFinderAlg
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
SmallClusterFinder(fhicl::ParameterSet const &pset)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
void FindSmallClusters(util::GeometryUtilities const &gser, detinfo::DetectorClocksData const &dataClocks, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> allHits)
Cluster finding and building.
static const SentryArgument_t Sentry
An instance of the sentry object.
Definition: Cluster.h:174
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:155
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
geo::WireReadoutGeom const & wireReadoutGeom
Interface for a class providing readout channel mapping to geometry.
Helper functions to create a cluster.
Wrapper for ClusterParamsAlgBase objects to accept diverse input.
std::vector< art::Ptr< recob::Hit > > GetLeftoversByPlane(unsigned int iPlane)
Wrapper for ClusterParamsAlgBase objects to accept arbitrary input.
Declaration of cluster object.
unsigned int Nplanes(TPCID const &tpcid=details::tpc_zero) const
Returns the total number of planes in the specified TPC.
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.
int GetPlaneAndTPC(art::Ptr< recob::Hit > a, unsigned int &p, unsigned int &cs, unsigned int &t, unsigned int &w)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
Encapsulate the construction of a single detector plane .
std::vector< std::vector< art::Ptr< recob::Hit > > > GetSmallClustersByPlane(unsigned int iPlane)
art::ServiceHandle< geo::Geometry const > geom
Interface to class computing cluster parameters.
PlaneGeo const & Plane(TPCID const &tpcid, View_t view) const
Returns the specified wire.
TCEvent evt
Definition: DataStructs.cxx:8
Float_t w
Definition: plot.C:20
art framework interface to geometry description