LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DBCluster3D_module.cc
Go to the documentation of this file.
1 // Class: DBCluster3D
3 // Plugin Type: producer (art v2_11_02)
4 // File: DBCluster3D_module.cc
5 //
6 // Generated at Sat Jun 16 07:43:32 2018 by Tingjun Yang using cetskelgen
7 // from cetlib version v3_03_01.
9 
16 #include "fhiclcpp/ParameterSet.h"
17 
26 
27 #include <memory>
28 
29 namespace cluster {
30  class DBCluster3D;
31 }
32 
34 
37 
38 public:
39  explicit DBCluster3D(fhicl::ParameterSet const& p);
40  // The compiler-generated destructor is fine for non-base
41  // classes without bare pointers or other resource use.
42 
43  // Plugins should not be copied or assigned.
44  DBCluster3D(DBCluster3D const&) = delete;
45  DBCluster3D(DBCluster3D&&) = delete;
46  DBCluster3D& operator=(DBCluster3D const&) = delete;
47  DBCluster3D& operator=(DBCluster3D&&) = delete;
48 
49 private:
50  // Required functions.
51  void produce(art::Event& e) override;
52 
56 
58 
60 
61  double tickToDist;
62  double fMinHitDis;
63 };
64 
66  : EDProducer{p}
67  , fHitModuleLabel(p.get<art::InputTag>("HitModuleLabel"))
68  , fSpacePointModuleLabel(p.get<art::InputTag>("SpacePointModuleLabel"))
69  , fSPHitAssnLabel(p.get<art::InputTag>("SPHitAssnLabel"))
70  , fDBScan(p.get<fhicl::ParameterSet>("DBScan3DAlg"))
71  , fMinHitDis(p.get<double>("MinHitDis"))
72 {
73  produces<std::vector<recob::Slice>>();
74  produces<art::Assns<recob::Slice, recob::Hit>>();
75  produces<art::Assns<recob::Slice, recob::SpacePoint>>();
76 
78  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
79  auto const det_prop =
81 
82  tickToDist = det_prop.DriftVelocity(det_prop.Efield(), det_prop.Temperature());
83  tickToDist *= 1.e-3 * sampling_rate(clock_data); // 1e-3 is conversion of 1/us to 1/ns
85 }
86 
88 {
89  auto scol = std::make_unique<std::vector<recob::Slice>>();
90  auto slc_hit_assn = std::make_unique<art::Assns<recob::Slice, recob::Hit>>();
91  auto slc_sps_assn = std::make_unique<art::Assns<recob::Slice, recob::SpacePoint>>();
92 
93  auto hitsHandle = evt.getValidHandle<std::vector<recob::Hit>>(fHitModuleLabel);
94  auto spsHandle = evt.getValidHandle<std::vector<recob::SpacePoint>>(fSpacePointModuleLabel);
95 
96  // all hits in the collection
97  std::vector<art::Ptr<recob::Hit>> hits;
98  art::fill_ptr_vector(hits, hitsHandle);
99 
100  // all space points in the collection
101  std::vector<art::Ptr<recob::SpacePoint>> sps;
102  art::fill_ptr_vector(sps, spsHandle);
103 
104  art::FindManyP<recob::SpacePoint> spFromHit(hitsHandle, evt, fSPHitAssnLabel);
105  if (!spFromHit.isValid()) {
106  std::cout << "spFromHit is invalid\n";
107  return;
108  }
109  // Find the first Hit - SpacePoint assn and check consistency on the first event
110  static bool first = true;
111  if (first) {
112  bool success = false;
113  bool foundsps = false;
114  for (auto& hit : hits) {
115  auto& sps = spFromHit.at(hit.key());
116  if (sps.empty()) continue;
117  success = (sps[0].id() == spsHandle.id());
118  foundsps = true;
119  break;
120  } // hit
121  if ((!success) && foundsps)
122  throw cet::exception("DBCluster3D")
123  << "HitModuleLabel, SpacePointModuleLabel and SPHitAssnLabel are inconsistent\n";
124  first = false;
125  } // first
126 
127  art::FindManyP<recob::Hit> hitFromSp(spsHandle, evt, fSPHitAssnLabel);
128  if (!hitFromSp.isValid()) {
129  std::cout << "hitFromSp is invalid\n";
130  return;
131  }
132  fDBScan.init(sps, hitFromSp);
133  fDBScan.dbscan();
134 
135  //Find number of slices
136  int maxid = INT_MIN;
137  for (size_t i = 0; i < fDBScan.points.size(); ++i) {
138  if (fDBScan.points[i].cluster_id > maxid) maxid = fDBScan.points[i].cluster_id;
139  }
140  size_t nslc = 0;
141  if (maxid >= 0) nslc = maxid + 1;
142 
143  //Save hits associated with each slice
144  std::vector<std::vector<art::Ptr<recob::Hit>>> slcHits(nslc);
145  //Save hits on each PlaneID with pfparticle index
146  std::map<geo::PlaneID, std::vector<std::pair<art::Ptr<recob::Hit>, unsigned int>>> hitmap;
147  for (auto& hit : hits) {
148  auto& sps = spFromHit.at(hit.key());
149  if (sps.size()) { //found associated space point
150  if (fDBScan.points[sps[0].key()].cluster_id >= 0) {
151  slcHits[fDBScan.points[sps[0].key()].cluster_id].push_back(hit);
152  hitmap[geo::PlaneID(hit->WireID())].push_back(
153  std::make_pair(hit, fDBScan.points[sps[0].key()].cluster_id));
154  }
155  } // sps.size()
156  } // hit
157 
158  //Save hits not associated with any spacepoints
159  for (auto& hit : hits) {
160  bool found = false;
161  for (size_t i = 0; i < slcHits.size(); ++i) {
162  if (std::find(slcHits[i].begin(), slcHits[i].end(), hit) != slcHits[i].end()) {
163  found = true;
164  break;
165  }
166  }
167  if (!found) {
168  double wirePitch = fGeom->WirePitch(hit->WireID());
169  double UnitsPerTick = tickToDist / wirePitch;
170  double x0 = hit->WireID().Wire;
171  double y0 = hit->PeakTime() * UnitsPerTick;
172  double mindis = DBL_MAX;
173  unsigned slcIndex = UINT_MAX;
174  for (auto& hit2 : hitmap[geo::PlaneID(hit->WireID())]) {
175  double dx = hit2.first->WireID().Wire - x0;
176  double dy = hit2.first->PeakTime() * UnitsPerTick - y0;
177  double dis = dx * dx + dy * dy;
178  if (dis < mindis) {
179  mindis = dis;
180  slcIndex = hit2.second;
181  }
182  }
183  if (slcIndex != UINT_MAX && mindis < fMinHitDis) { slcHits[slcIndex].push_back(hit); }
184  }
185  }
186 
187  //Save spacepoints for each slice
188  std::vector<std::vector<art::Ptr<recob::SpacePoint>>> sps_in_slc(nslc);
189  for (size_t i = 0; i < fDBScan.points.size(); ++i) {
190  if (fDBScan.points[i].cluster_id >= 0) {
191  sps_in_slc[fDBScan.points[i].cluster_id].push_back(sps[i]);
192  }
193  } // i
194 
195  // calculate the properties of the slice
196  for (size_t isl = 0; isl < nslc; ++isl) {
197  double sum = sps_in_slc[isl].size();
198  // find the center
199  double center[3] = {0.};
200  // TODO: calculate the charge using recob::PointCharge instead of just counting hits
201  float charge = slcHits[isl].size();
202  for (auto& spt : sps_in_slc[isl]) {
203  for (unsigned short xyz = 0; xyz < 3; ++xyz)
204  center[xyz] += spt->XYZ()[xyz];
205  } // spt
206  for (unsigned short xyz = 0; xyz < 3; ++xyz)
207  center[xyz] /= sum;
208  // do a linear fit
209  double sumx = 0, sumy = 0., sumz = 0., sumx2 = 0, sumy2 = 0, sumz2 = 0;
210  double sumxy = 0, sumxz = 0;
211  for (auto& spt : sps_in_slc[isl]) {
212  double xx = spt->XYZ()[0] - center[0];
213  double yy = spt->XYZ()[1] - center[1];
214  double zz = spt->XYZ()[2] - center[2];
215  sumx += xx;
216  sumy += yy;
217  sumz += zz;
218  sumx2 += xx * xx;
219  sumy2 += yy * yy;
220  sumz2 += zz * zz;
221  sumxy += xx * yy;
222  sumxz += xx * zz;
223  } // spt
224  double delta = sum * sumx2 - sumx * sumx;
225  if (delta <= 0) continue;
226  // calculate the slopes
227  double dydx = (sumxy * sum - sumx * sumy) / delta;
228  double dzdx = (sumxz * sum - sumx * sumz) / delta;
229  // and convert to direction vector
230  double direction[3];
231  double norm = std::sqrt(1 + dydx * dydx + dzdx * dzdx);
232  direction[0] = 1 / norm;
233  direction[1] = dydx / norm;
234  direction[2] = dzdx / norm;
235  // Find the points that are farthest away from the center along the slice axis
236  unsigned int imax = 0, imin = 0;
237  float minAlong = 1E6, maxAlong = -1E6;
238  double tmpVec[3];
239  for (unsigned int ipt = 0; ipt < sps_in_slc[isl].size(); ++ipt) {
240  auto& spt = sps_in_slc[isl][ipt];
241  for (unsigned short xyz = 0; xyz < 3; ++xyz)
242  tmpVec[xyz] = spt->XYZ()[xyz] - center[xyz];
243  double dotp = 0;
244  for (unsigned short xyz = 0; xyz < 3; ++xyz)
245  dotp += tmpVec[xyz] * direction[xyz];
246  if (dotp < minAlong) {
247  minAlong = dotp;
248  imin = ipt;
249  }
250  if (dotp > maxAlong) {
251  maxAlong = dotp;
252  imax = ipt;
253  }
254  } // spt
255  // Find the aspect ratio
256  float cnt = 0.;
257  float aspectRatio = 0;
258  double arg = sum * sumy2 - sumy * sumy;
259  if (arg > 0) {
260  aspectRatio += std::abs(sum * sumxy - sumx * sumy) / sqrt(delta * arg);
261  ++cnt;
262  }
263  arg = sum * sumz2 - sumz * sumz;
264  if (arg > 0) {
265  aspectRatio += std::abs(sum * sumxz - sumx * sumz) / sqrt(delta * arg);
266  ++cnt;
267  }
268  if (cnt > 1) aspectRatio /= cnt;
269  int id = isl + 1;
270  Point_t ctr(center[0], center[1], center[2]);
271  Vector_t dir(direction[0], direction[1], direction[2]);
272  auto pos0 = sps_in_slc[isl][imin]->XYZ();
273  Point_t ep0(pos0[0], pos0[1], pos0[2]);
274  auto pos1 = sps_in_slc[isl][imax]->XYZ();
275  Point_t ep1(pos1[0], pos1[1], pos1[2]);
276  scol->emplace_back(id, ctr, dir, ep0, ep1, aspectRatio, charge);
277  util::CreateAssn(evt, *scol, slcHits[isl], *slc_hit_assn);
278  util::CreateAssn(evt, *scol, sps_in_slc[isl], *slc_sps_assn);
279  } // isl
280 
281  evt.put(std::move(scol));
282  evt.put(std::move(slc_hit_assn));
283  evt.put(std::move(slc_sps_assn));
284 }
285 
geo::GeometryCore const * fGeom
DBCluster3D & operator=(DBCluster3D const &)=delete
Double_t xx
Definition: macro.C:12
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:463
constexpr auto abs(T v)
Returns the absolute value of the argument.
Double_t zz
Definition: plot.C:276
const art::InputTag fSpacePointModuleLabel
Cluster finding and building.
const art::InputTag fHitModuleLabel
const art::InputTag fSPHitAssnLabel
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
recob::tracking::Vector_t Vector_t
void hits()
Definition: readHits.C:15
DBCluster3D(fhicl::ParameterSet const &p)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void produce(art::Event &e) override
void init(const std::vector< art::Ptr< recob::SpacePoint >> &sps, art::FindManyP< recob::Hit > &hitFromSp)
Definition: DBScan3DAlg.cxx:29
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space. See recob::tracking::Coord_t for more details on the ...
Definition: TrackingTypes.h:31
Description of geometry of one entire detector.
Definition: GeometryCore.h:119
recob::tracking::Point_t Point_t
Detector simulation of raw signals on wires.
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.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Utility object to perform functions of association.
Float_t norm
TDirectory * dir
Definition: macro.C:5
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
std::vector< point_t > points
Definition: DBScan3DAlg.h:82
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space. See recob::tracking::Coord_t for more detai...
Definition: TrackingTypes.h:27
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33