Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
cluster::DBCluster3D Class Reference
Inheritance diagram for cluster::DBCluster3D:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 DBCluster3D (fhicl::ParameterSet const &p)
 
 DBCluster3D (DBCluster3D const &)=delete
 
 DBCluster3D (DBCluster3D &&)=delete
 
DBCluster3D & operator= (DBCluster3D const &)=delete
 
DBCluster3D & operator= (DBCluster3D &&)=delete
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Producer & operator= (Producer const &)=delete
 
Producer & operator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifier & operator= (Modifier const &)=delete
 
Modifier & operator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Types

using Point_t = recob::tracking::Point_t
 
using Vector_t = recob::tracking::Vector_t
 

Private Member Functions

void produce (art::Event &e) override
 

Private Attributes

const art::InputTag fHitModuleLabel
 
const art::InputTag fSpacePointModuleLabel
 
const art::InputTag fSPHitAssnLabel
 
DBScan3DAlg fDBScan
 
geo::GeometryCore const * fGeom
 
double tickToDist
 
double fMinHitDis
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 32 of file DBCluster3D_module.cc.

Member Typedef Documentation

using cluster::DBCluster3D::Point_t = recob::tracking::Point_t
private

Definition at line 34 of file DBCluster3D_module.cc.

using cluster::DBCluster3D::Vector_t = recob::tracking::Vector_t
private

Definition at line 35 of file DBCluster3D_module.cc.

Constructor & Destructor Documentation

cluster::DBCluster3D::DBCluster3D ( fhicl::ParameterSet const &  p)
explicit

Definition at line 64 of file DBCluster3D_module.cc.

65  : EDProducer{p}
66  , fHitModuleLabel(p.get<art::InputTag>("HitModuleLabel"))
67  , fSpacePointModuleLabel(p.get<art::InputTag>("SpacePointModuleLabel"))
68  , fSPHitAssnLabel(p.get<art::InputTag>("SPHitAssnLabel"))
69  , fDBScan(p.get<fhicl::ParameterSet>("DBScan3DAlg"))
70  , fMinHitDis(p.get<double>("MinHitDis"))
71 {
72  produces<std::vector<recob::Slice>>();
73  produces<art::Assns<recob::Slice, recob::Hit>>();
74  produces<art::Assns<recob::Slice, recob::SpacePoint>>();
75 
76  fGeom = art::ServiceHandle<geo::Geometry const>().get();
77  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
78  auto const det_prop =
79  art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clock_data);
80 
81  tickToDist = det_prop.DriftVelocity(det_prop.Efield(), det_prop.Temperature());
82  tickToDist *= 1.e-3 * sampling_rate(clock_data); // 1e-3 is conversion of 1/us to 1/ns
83  fMinHitDis *= fMinHitDis;
84 }
cluster::DBCluster3D::fGeom
geo::GeometryCore const * fGeom
Definition: DBCluster3D_module.cc:58
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
cluster::DBCluster3D::fSpacePointModuleLabel
const art::InputTag fSpacePointModuleLabel
Definition: DBCluster3D_module.cc:53
cluster::DBCluster3D::fHitModuleLabel
const art::InputTag fHitModuleLabel
Definition: DBCluster3D_module.cc:52
cluster::DBCluster3D::fSPHitAssnLabel
const art::InputTag fSPHitAssnLabel
Definition: DBCluster3D_module.cc:54
test.p
p
Definition: test.py:223
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
cluster::DBCluster3D::DBCluster3D ( DBCluster3D const &  )
delete
cluster::DBCluster3D::DBCluster3D ( DBCluster3D &&  )
delete

Member Function Documentation

DBCluster3D& cluster::DBCluster3D::operator= ( DBCluster3D const &  )
delete
DBCluster3D& cluster::DBCluster3D::operator= ( DBCluster3D &&  )
delete
void cluster::DBCluster3D::produce ( art::Event &  e)
overrideprivatevirtual

Implements art::EDProducer.

Definition at line 87 of file DBCluster3D_module.cc.

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 }
cluster::DBCluster3D::fGeom
geo::GeometryCore const * fGeom
Definition: DBCluster3D_module.cc:58
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
Point_t
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >> Point_t
Definition: EnergyRecoOutput.h:8
geo::PlaneID
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
cluster::DBCluster3D::fSpacePointModuleLabel
const art::InputTag fSpacePointModuleLabel
Definition: DBCluster3D_module.cc:53
dir
string dir
Definition: INukeNucleonCorr.cxx:61
cluster::DBCluster3D::fHitModuleLabel
const art::InputTag fHitModuleLabel
Definition: DBCluster3D_module.cc:52
cluster::DBCluster3D::fSPHitAssnLabel
const art::InputTag fSPHitAssnLabel
Definition: DBCluster3D_module.cc:54
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:824
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
push_back
fInnerVessel push_back(Point(-578.400000, 0.000000, 0.000000))
cluster::DBScan3DAlg::init
void init(const std::vector< art::Ptr< recob::SpacePoint >> &sps, art::FindManyP< recob::Hit > &hitFromSp)
Definition: DBScan3DAlg.cxx:29
wirecell.gen.depos.move
def move(depos, offset)
Definition: depos.py:107
util::CreateAssn
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
Definition: AssociationUtil.h:601
geo::vect::norm
auto norm(Vector const &v)
Return norm of the specified vector.
Definition: geo_vectors_utils.h:997
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
tracks.imax
int imax
Definition: tracks.py:195
wirecell.gen.depos.center
def center(depos, point)
Definition: depos.py:117
CRT::PlaneID
detail::Node< FrameID, bool > PlaneID
Definition: CRTID.h:125
util::begin
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
cluster::DBScan3DAlg::points
std::vector< point_t > points
Definition: DBScan3DAlg.h:83
tca::evt
TCEvent evt
Definition: DataStructs.cxx:7
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
Vector_t
recob::tracking::Vector_t Vector_t
Definition: PMAlgTracking.cxx:21
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

DBScan3DAlg cluster::DBCluster3D::fDBScan
private

Definition at line 56 of file DBCluster3D_module.cc.

geo::GeometryCore const* cluster::DBCluster3D::fGeom
private

Definition at line 58 of file DBCluster3D_module.cc.

const art::InputTag cluster::DBCluster3D::fHitModuleLabel
private

Definition at line 52 of file DBCluster3D_module.cc.

double cluster::DBCluster3D::fMinHitDis
private

Definition at line 61 of file DBCluster3D_module.cc.

const art::InputTag cluster::DBCluster3D::fSpacePointModuleLabel
private

Definition at line 53 of file DBCluster3D_module.cc.

const art::InputTag cluster::DBCluster3D::fSPHitAssnLabel
private

Definition at line 54 of file DBCluster3D_module.cc.

double cluster::DBCluster3D::tickToDist
private

Definition at line 60 of file DBCluster3D_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.11