LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DBcluster_module.cc
Go to the documentation of this file.
1 //
3 // \file DBcluster_module.cc
4 //
5 // \author kinga.partyka@yale.edu
6 //
8 
9 //Framework includes:
15 #include "art_root_io/TFileService.h"
18 #include "fhiclcpp/ParameterSet.h"
20 
21 //LArSoft includes
29 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
30 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
35 
36 #include "TH1.h"
37 #include <cstdlib>
38 #include <iomanip>
39 
40 namespace cluster {
41 
42  //---------------------------------------------------------------
43  class DBcluster : public art::EDProducer {
44  public:
45  explicit DBcluster(fhicl::ParameterSet const& pset);
46 
47  private:
48  void produce(art::Event& evt);
49  void beginJob();
50 
51  TH1F* fhitwidth;
54 
55  std::string fhitsModuleLabel;
56 
58  };
59 
60 }
61 
62 namespace cluster {
63 
64  //-------------------------------------------------
66  : EDProducer{pset}, fDBScan(pset.get<fhicl::ParameterSet>("DBScanAlg"))
67  {
68  fhitsModuleLabel = pset.get<std::string>("HitsModuleLabel");
69 
70  produces<std::vector<recob::Cluster>>();
71  produces<art::Assns<recob::Cluster, recob::Hit>>();
72  }
73 
74  //-------------------------------------------------
76  {
77  // get access to the TFile service
79 
80  fhitwidth = tfs->make<TH1F>(" fhitwidth", "width of hits in cm", 50000, 0, 5);
81  fhitwidth_ind_test = tfs->make<TH1F>("fhitwidth_ind_test", "width of hits in cm", 50000, 0, 5);
83  tfs->make<TH1F>("fhitwidth_coll_test", "width of hits in cm", 50000, 0, 5);
84  }
85 
86  //-----------------------------------------------------------------
88  {
89 
90  //get a collection of clusters
91  std::unique_ptr<std::vector<recob::Cluster>> ccol(new std::vector<recob::Cluster>);
92  std::unique_ptr<art::Assns<recob::Cluster, recob::Hit>> assn(
94 
95  // prepare the algorithm to compute the cluster characteristics;
96  // we use the "standard" one here; configuration would happen here,
97  // but we are using the default configuration for that algorithm
99 
101 
103  evt.getByLabel(fhitsModuleLabel, hitcol);
104 
105  // loop over all hits in the event and look for clusters (for each plane)
106  std::vector<art::Ptr<recob::Hit>> allhits;
107 
108  // get channel quality service:
109  lariov::ChannelStatusProvider const& channelStatus =
111 
112  lariov::ChannelStatusProvider::ChannelSet_t const BadChannels = channelStatus.BadChannels();
113 
114  // make a map of the geo::PlaneID to vectors of art::Ptr<recob::Hit>
115  std::map<geo::PlaneID, std::vector<art::Ptr<recob::Hit>>> planeIDToHits;
116  for (size_t i = 0; i < hitcol->size(); ++i)
117  planeIDToHits[hitcol->at(i).WireID().planeID()].push_back(art::Ptr<recob::Hit>(hitcol, i));
118 
119  auto const clock_data =
121  auto const det_prop =
123  util::GeometryUtilities const gser{*geom, clock_data, det_prop};
124  for (auto& itr : planeIDToHits) {
125 
126  geo::SigType_t sigType = geom->SignalType(itr.first);
127  allhits.resize(itr.second.size());
128  allhits.swap(itr.second);
129 
130  fDBScan.InitScan(clock_data, det_prop, allhits, BadChannels);
131 
132  //----------------------------------------------------------------
133  for (unsigned int j = 0; j < fDBScan.fps.size(); ++j) {
134 
135  if (allhits.size() != fDBScan.fps.size()) break;
136 
137  fhitwidth->Fill(fDBScan.fps[j][2]);
138 
139  if (sigType == geo::kInduction) fhitwidth_ind_test->Fill(fDBScan.fps[j][2]);
140  if (sigType == geo::kCollection) fhitwidth_coll_test->Fill(fDBScan.fps[j][2]);
141  }
142 
143  //*******************************************************************
145 
146  for (size_t i = 0; i < fDBScan.fclusters.size(); ++i) {
147  art::PtrVector<recob::Hit> clusterHits;
148 
149  for (size_t j = 0; j < fDBScan.fpointId_to_clusterId.size(); ++j) {
150  if (fDBScan.fpointId_to_clusterId[j] == i) { clusterHits.push_back(allhits[j]); }
151  }
152 
153  if (clusterHits.empty()) continue;
154 
156  const geo::WireID& wireID = clusterHits.front()->WireID();
157  unsigned int sw = wireID.Wire;
158  unsigned int ew = clusterHits.back()->WireID().Wire;
159 
160  // feed the algorithm with all the cluster hits
161  ClusterParamAlgo.ImportHits(gser, clusterHits);
162 
163  // create the recob::Cluster directly in the vector
164  ClusterCreator cluster(gser,
165  ClusterParamAlgo, // algo
166  float(sw), // start_wire
167  0., // sigma_start_wire
168  clusterHits.front()->PeakTime(), // start_tick
169  clusterHits.front()->SigmaPeakTime(), // sigma_start_tick
170  float(ew), // end_wire
171  0., // sigma_end_wire,
172  clusterHits.back()->PeakTime(), // end_tick
173  clusterHits.back()->SigmaPeakTime(), // sigma_end_tick
174  ccol->size(), // ID
175  clusterHits.front()->View(), // view
176  wireID.planeID(), // plane
177  recob::Cluster::Sentry // sentry
178  );
179 
180  ccol->emplace_back(cluster.move());
181 
182  // associate the hits to this cluster
183  util::CreateAssn(evt, *ccol, clusterHits, *assn);
184 
185  } //end loop over fclusters
186 
187  allhits.clear();
188  } // end loop over PlaneID map
189 
190  mf::LogVerbatim("Summary") << std::setfill('-') << std::setw(175) << "-" << std::setfill(' ');
191  mf::LogVerbatim("Summary") << "DBcluster Summary:";
192  for (unsigned int i = 0; i < ccol->size(); ++i)
193  mf::LogVerbatim("Summary") << ccol->at(i);
194 
195  evt.put(std::move(ccol));
196  evt.put(std::move(assn));
197 
198  return;
199  }
200 
201 } // end namespace
202 
203 namespace cluster {
204 
206 
207 }
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Class managing the creation of a new recob::Cluster object.
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
std::vector< std::vector< double > > fps
the collection of points we are working on
Definition: DBScanAlg.h:71
Cluster finding and building.
static const SentryArgument_t Sentry
An instance of the sentry object.
Definition: Cluster.h:174
DBScanAlg fDBScan
object that implements the DB scan algorithm
void ImportHits(util::GeometryUtilities const &gser, Iter begin, Iter end)
Calls SetHits() with the hits in the sequence.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
reference back()
Definition: PtrVector.h:387
std::vector< unsigned int > fpointId_to_clusterId
mapping point_id -> clusterId
Definition: DBScanAlg.h:72
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
Signal from induction planes.
Definition: geo_types.h:151
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
Helper functions to create a cluster.
Wrapper for ClusterParamsAlgBase objects to accept diverse input.
Wrapper for ClusterParamsAlgBase objects to accept arbitrary input.
Declaration of cluster object.
bool empty() const
Definition: PtrVector.h:330
std::string fhitsModuleLabel
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.
reference front()
Definition: PtrVector.h:373
void InitScan(const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &allhits, std::set< uint32_t > badChannels, const std::vector< geo::WireID > &wireids=std::vector< geo::WireID >())
Definition: DBScanAlg.cxx:255
Float_t sw
Definition: plot.C:20
void produce(art::Event &evt)
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
std::vector< std::vector< unsigned int > > fclusters
collection of something
Definition: DBScanAlg.h:70
DBcluster(fhicl::ParameterSet const &pset)
constexpr PlaneID const & planeID() const
Definition: geo_types.h:620
Interface to class computing cluster parameters.
TCEvent evt
Definition: DataStructs.cxx:8
art framework interface to geometry description
Signal from collection planes.
Definition: geo_types.h:152