LArSoft  v10_04_05
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
30 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
31 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
36 
37 #include "TH1.h"
38 #include <cstdlib>
39 #include <iomanip>
40 
41 namespace cluster {
42 
43  //---------------------------------------------------------------
44  class DBcluster : public art::EDProducer {
45  public:
46  explicit DBcluster(fhicl::ParameterSet const& pset);
47 
48  private:
49  void produce(art::Event& evt);
50  void beginJob();
51 
52  TH1F* fhitwidth;
55 
56  std::string fhitsModuleLabel;
57 
59  };
60 
61 }
62 
63 namespace cluster {
64 
65  //-------------------------------------------------
67  : EDProducer{pset}, fDBScan(pset.get<fhicl::ParameterSet>("DBScanAlg"))
68  {
69  fhitsModuleLabel = pset.get<std::string>("HitsModuleLabel");
70 
71  produces<std::vector<recob::Cluster>>();
72  produces<art::Assns<recob::Cluster, recob::Hit>>();
73  }
74 
75  //-------------------------------------------------
77  {
78  // get access to the TFile service
80 
81  fhitwidth = tfs->make<TH1F>(" fhitwidth", "width of hits in cm", 50000, 0, 5);
82  fhitwidth_ind_test = tfs->make<TH1F>("fhitwidth_ind_test", "width of hits in cm", 50000, 0, 5);
84  tfs->make<TH1F>("fhitwidth_coll_test", "width of hits in cm", 50000, 0, 5);
85  }
86 
87  //-----------------------------------------------------------------
89  {
90 
91  //get a collection of clusters
92  std::unique_ptr<std::vector<recob::Cluster>> ccol(new std::vector<recob::Cluster>);
93  std::unique_ptr<art::Assns<recob::Cluster, recob::Hit>> assn(
95 
96  // prepare the algorithm to compute the cluster characteristics;
97  // we use the "standard" one here; configuration would happen here,
98  // but we are using the default configuration for that algorithm
100 
102  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
103 
105  evt.getByLabel(fhitsModuleLabel, hitcol);
106 
107  // loop over all hits in the event and look for clusters (for each plane)
108  std::vector<art::Ptr<recob::Hit>> allhits;
109 
110  // get channel quality service:
111  lariov::ChannelStatusProvider const& channelStatus =
113 
114  lariov::ChannelStatusProvider::ChannelSet_t const BadChannels = channelStatus.BadChannels();
115 
116  // make a map of the geo::PlaneID to vectors of art::Ptr<recob::Hit>
117  std::map<geo::PlaneID, std::vector<art::Ptr<recob::Hit>>> planeIDToHits;
118  for (size_t i = 0; i < hitcol->size(); ++i)
119  planeIDToHits[hitcol->at(i).WireID().planeID()].push_back(art::Ptr<recob::Hit>(hitcol, i));
120 
121  auto const clock_data =
123  auto const det_prop =
125  util::GeometryUtilities const gser{*geom, wireReadoutGeom, clock_data, det_prop};
126  for (auto& itr : planeIDToHits) {
127 
128  geo::SigType_t sigType = wireReadoutGeom.SignalType(itr.first);
129  allhits.resize(itr.second.size());
130  allhits.swap(itr.second);
131 
132  fDBScan.InitScan(clock_data, det_prop, allhits, BadChannels);
133 
134  //----------------------------------------------------------------
135  for (unsigned int j = 0; j < fDBScan.fps.size(); ++j) {
136 
137  if (allhits.size() != fDBScan.fps.size()) break;
138 
139  fhitwidth->Fill(fDBScan.fps[j][2]);
140 
141  if (sigType == geo::kInduction) fhitwidth_ind_test->Fill(fDBScan.fps[j][2]);
142  if (sigType == geo::kCollection) fhitwidth_coll_test->Fill(fDBScan.fps[j][2]);
143  }
144 
145  //*******************************************************************
147 
148  for (size_t i = 0; i < fDBScan.fclusters.size(); ++i) {
149  art::PtrVector<recob::Hit> clusterHits;
150 
151  for (size_t j = 0; j < fDBScan.fpointId_to_clusterId.size(); ++j) {
152  if (fDBScan.fpointId_to_clusterId[j] == i) { clusterHits.push_back(allhits[j]); }
153  }
154 
155  if (clusterHits.empty()) continue;
156 
158  const geo::WireID& wireID = clusterHits.front()->WireID();
159  unsigned int sw = wireID.Wire;
160  unsigned int ew = clusterHits.back()->WireID().Wire;
161 
162  // feed the algorithm with all the cluster hits
163  ClusterParamAlgo.ImportHits(gser, clusterHits);
164 
165  // create the recob::Cluster directly in the vector
166  ClusterCreator cluster(gser,
167  ClusterParamAlgo, // algo
168  float(sw), // start_wire
169  0., // sigma_start_wire
170  clusterHits.front()->PeakTime(), // start_tick
171  clusterHits.front()->SigmaPeakTime(), // sigma_start_tick
172  float(ew), // end_wire
173  0., // sigma_end_wire,
174  clusterHits.back()->PeakTime(), // end_tick
175  clusterHits.back()->SigmaPeakTime(), // sigma_end_tick
176  ccol->size(), // ID
177  clusterHits.front()->View(), // view
178  wireID.planeID(), // plane
179  recob::Cluster::Sentry // sentry
180  );
181 
182  ccol->emplace_back(cluster.move());
183 
184  // associate the hits to this cluster
185  util::CreateAssn(evt, *ccol, clusterHits, *assn);
186 
187  } //end loop over fclusters
188 
189  allhits.clear();
190  } // end loop over PlaneID map
191 
192  mf::LogVerbatim("Summary") << std::setfill('-') << std::setw(175) << "-" << std::setfill(' ');
193  mf::LogVerbatim("Summary") << "DBcluster Summary:";
194  for (unsigned int i = 0; i < ccol->size(); ++i)
195  mf::LogVerbatim("Summary") << ccol->at(i);
196 
197  evt.put(std::move(ccol));
198  evt.put(std::move(assn));
199  }
200 
201 } // end namespace
202 
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
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
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:147
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:256
Float_t sw
Definition: plot.C:20
void produce(art::Event &evt)
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:471
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:148