VertexFinder2D_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // VertexFinder2D class
4 //
5 // tjyang@fnal.gov
6 //
7 // This algorithm is designed to reconstruct the vertices using the
8 // 2D cluster information
9 //
10 // This is Preliminary Work and needs modifications
11 // ////////////////////////////////////////////////////////////////////////
12 #include <algorithm>
13 #include <cmath>
14 #include <iomanip>
15 #include <string>
16 
17 // Framework includes
23 #include "art_root_io/TFileService.h"
24 #include "canvas/Persistency/Common/FindManyP.h"
27 #include "fhiclcpp/ParameterSet.h"
29 
41 
42 #include "TF1.h"
43 #include "TGraph.h"
44 #include "TH1D.h"
45 #include "TMath.h"
46 
47 namespace {
48  struct CluLen {
49  int index;
50  float length;
51  };
52 
53  bool
54  myfunction(CluLen c1, CluLen c2)
55  {
56  return (c1.length > c2.length);
57  }
58 
59  struct SortByWire {
60  bool
61  operator()(art::Ptr<recob::Hit> const& h1, art::Ptr<recob::Hit> const& h2) const
62  {
63  return h1->Channel() < h2->Channel();
64  }
65  };
66 }
67 
68 ///vertex reconstruction
69 namespace vertex {
70 
72  public:
73  explicit VertexFinder2D(fhicl::ParameterSet const& pset);
74 
75  private:
76  void beginJob() override;
77  void produce(art::Event& evt) override;
78 
79  TH1D* dtIC;
80 
82  };
83 
84 }
85 
86 namespace vertex {
87 
88  //-----------------------------------------------------------------------------
89  VertexFinder2D::VertexFinder2D(fhicl::ParameterSet const& pset) : EDProducer{pset}
90  {
91  fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
92  produces<std::vector<recob::Vertex>>();
93  produces<std::vector<recob::EndPoint2D>>();
94  produces<art::Assns<recob::EndPoint2D, recob::Hit>>();
95  produces<art::Assns<recob::Vertex, recob::Hit>>();
96  produces<art::Assns<recob::Vertex, recob::Shower>>();
97  produces<art::Assns<recob::Vertex, recob::Track>>();
98  }
99 
100  //-------------------------------------------------------------------------
101  void
103  {
104  // get access to the TFile service
106  dtIC = tfs->make<TH1D>("dtIC", "It0-Ct0", 100, -5, 5);
107  dtIC->Sumw2();
108  }
109 
110  // //-----------------------------------------------------------------------------
111  void
113  {
115  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
116  auto const detProp =
118 
119  // define TPC parameters
120  TString tpcName = geom->GetLArTPCVolumeName();
121 
122  double YC = (geom->DetHalfHeight()) * 2.;
123 
124  // wire angle with respect to the vertical direction
125  double Angle = geom->Plane(1).Wire(0).ThetaZ(false) - TMath::Pi() / 2.;
126 
127  // Parameters temporary defined here, but possibly to be retrieved somewhere
128  // in the code
129  double timetick = sampling_rate(clockData) * 1.e-3; // time sample in us
130  double presamplings = trigger_offset(clockData);
131 
132  double wire_pitch = geom->WirePitch(); //wire pitch in cm
133  double Efield_drift = detProp.Efield(); // Electric Field in the drift region in kV/cm
134  double Temperature = detProp.Temperature(); // LAr Temperature in K
135 
136  //drift velocity in the drift region (cm/us)
137  double driftvelocity = detProp.DriftVelocity(Efield_drift, Temperature);
138 
139  //time sample (cm)
140  double timepitch = driftvelocity * timetick;
141 
142  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
143  evt.getByLabel(fClusterModuleLabel, clusterListHandle);
144 
146  for (unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
147  art::Ptr<recob::Cluster> clusterHolder(clusterListHandle, ii);
148  clusters.push_back(clusterHolder);
149  }
150 
151  art::FindManyP<recob::Hit> fmh(clusterListHandle, evt, fClusterModuleLabel);
152 
153  //Point to a collection of vertices to output.
154  std::unique_ptr<std::vector<recob::Vertex>> vcol(new std::vector<recob::Vertex>); //3D vertex
155  std::unique_ptr<std::vector<recob::EndPoint2D>> epcol(
156  new std::vector<recob::EndPoint2D>); //2D vertex
157  std::unique_ptr<art::Assns<recob::EndPoint2D, recob::Hit>> assnep(
159  std::unique_ptr<art::Assns<recob::Vertex, recob::Shower>> assnsh(
161  std::unique_ptr<art::Assns<recob::Vertex, recob::Track>> assntr(
163  std::unique_ptr<art::Assns<recob::Vertex, recob::Hit>> assnh(
165 
166  // nplanes here is really being used as a proxy for the
167  // number of views in the detector
168  int nplanes = geom->Views().size();
169 
170  std::vector<std::vector<int>> Cls(nplanes); //index to clusters in each view
171  std::vector<std::vector<CluLen>> clulens(nplanes);
172 
173  std::vector<double> dtdwstart;
174 
175  //loop over clusters
176  for (size_t iclu = 0; iclu < clusters.size(); ++iclu) {
177 
178  float w0 = clusters[iclu]->StartWire();
179  float w1 = clusters[iclu]->EndWire();
180  float t0 = clusters[iclu]->StartTick();
181  float t1 = clusters[iclu]->EndTick();
182 
183  CluLen clulen;
184  clulen.index = iclu;
185  clulen.length = std::hypot((w0 - w1) * wire_pitch,
186  detProp.ConvertTicksToX(t0, clusters[iclu]->View(), 0, 0) -
187  detProp.ConvertTicksToX(t1, clusters[iclu]->View(), 0, 0));
188 
189  switch (clusters[iclu]->View()) {
190 
191  case geo::kU: clulens[0].push_back(clulen); break;
192  case geo::kV: clulens[1].push_back(clulen); break;
193  case geo::kZ: clulens[2].push_back(clulen); break;
194  default: break;
195  }
196 
197  std::vector<double> wires;
198  std::vector<double> times;
199 
200  std::vector<art::Ptr<recob::Hit>> hit = fmh.at(iclu);
201  std::sort(hit.begin(), hit.end(), SortByWire());
202  int n = 0;
203  for (size_t i = 0; i < hit.size(); ++i) {
204  wires.push_back(hit[i]->WireID().Wire);
205  times.push_back(hit[i]->PeakTime());
206  ++n;
207  }
208  if (n >= 2) {
209  TGraph* the2Dtrack = new TGraph(std::min(10, n), &wires[0], &times[0]);
210  try {
211  the2Dtrack->Fit("pol1", "Q");
212  TF1* pol1 = (TF1*)the2Dtrack->GetFunction("pol1");
213  double par[2];
214  pol1->GetParameters(par);
215  dtdwstart.push_back(par[1]);
216  }
217  catch (...) {
218  mf::LogWarning("VertexFinder2D") << "Fitter failed";
219  delete the2Dtrack;
220  dtdwstart.push_back(std::tan(clusters[iclu]->StartAngle()));
221  continue;
222  }
223  delete the2Dtrack;
224  }
225  else
226  dtdwstart.push_back(std::tan(clusters[iclu]->StartAngle()));
227  }
228 
229  //sort clusters based on 2D length
230  for (size_t i = 0; i < clulens.size(); ++i) {
231  std::sort(clulens[i].begin(), clulens[i].end(), myfunction);
232  for (size_t j = 0; j < clulens[i].size(); ++j) {
233  Cls[i].push_back(clulens[i][j].index);
234  }
235  }
236 
237  std::vector<std::vector<int>> cluvtx(nplanes);
238  std::vector<double> vtx_w;
239  std::vector<double> vtx_t;
240 
241  for (int i = 0; i < nplanes; ++i) {
242  if (Cls[i].size() >= 1) {
243  //at least one cluster
244  //find the longest two clusters
245  int c1 = -1;
246  int c2 = -1;
247  double ww0 = -999;
248  double wb1 = -999;
249  double we1 = -999;
250  double wb2 = -999;
251  double we2 = -999;
252  double tt1 = -999;
253  double tt2 = -999;
254  double dtdw1 = -999;
255  double dtdw2 = -999;
256  double lclu1 = -999;
257  double lclu2 = -999;
258  for (unsigned j = 0; j < Cls[i].size(); ++j) {
259  double lclu = std::sqrt(
260  pow((clusters[Cls[i][j]]->StartWire() - clusters[Cls[i][j]]->EndWire()) * 13.5, 2) +
261  pow(clusters[Cls[i][j]]->StartTick() - clusters[Cls[i][j]]->EndTick(), 2));
262  bool rev = false;
263  bool deltaraylike = false;
264  bool enoughhits = false;
265  if (c1 != -1) {
266  double wb = clusters[Cls[i][j]]->StartWire();
267  double we = clusters[Cls[i][j]]->EndWire();
268  double tt = clusters[Cls[i][j]]->StartTick();
269  double dtdw = dtdwstart[Cls[i][j]];
270  int nhits = fmh.at(Cls[i][j]).size();
271  ww0 = (tt - tt1 + dtdw1 * wb1 - dtdw * wb) / (dtdw1 - dtdw);
272  if (std::abs(wb1 - ww0) > std::abs(we1 - ww0)) rev = true; //reverse cluster dir
273  if ((!rev && ww0 > wb1 + 15) || (rev && ww0 < we1 - 15)) deltaraylike = true;
274  if (((!rev && ww0 > wb1 + 10) || (rev && ww0 < we1 - 10)) && nhits < 5)
275  deltaraylike = true;
276  if (wb > wb1 + 20 && nhits < 20) deltaraylike = true;
277  if (wb > wb1 + 50 && nhits < 20) deltaraylike = true;
278  if (wb > wb1 + 8 && TMath::Abs(dtdw1 - dtdw) < 0.15) deltaraylike = true;
279  if (std::abs(wb - wb1) > 30 && std::abs(we - we1) > 30) deltaraylike = true;
280  if (std::abs(tt - tt1) > 100)
281  deltaraylike = true; //not really deltaray, but isolated cluster
282  //make sure there are enough hits in the cluster
283  //at leaset 2 hits if goes horizentally, at leaset 4 hits if goes vertically
284  double alpha = std::atan(dtdw);
285  if (nhits >= int(2 + 3 * (1 - std::abs(std::cos(alpha))))) enoughhits = true;
286  if (nhits < 5 && (ww0 < wb1 - 20 || ww0 > we1 + 20)) enoughhits = false;
287  }
288  //do not replace the second cluster if the 3rd cluster is not consistent with the existing 2
289  bool replace = true;
290  if (c1 != -1 && c2 != -1) {
291  double wb = clusters[Cls[i][j]]->StartWire();
292  double we = clusters[Cls[i][j]]->EndWire();
293  ww0 = (tt2 - tt1 + dtdw1 * wb1 - dtdw2 * wb2) / (dtdw1 - dtdw2);
294  if ((std::abs(ww0 - wb1) < 10 || std::abs(ww0 - we1) < 10) &&
295  (std::abs(ww0 - wb2) < 10 || std::abs(ww0 - we2) < 10)) {
296  if (std::abs(ww0 - wb) > 15 && std::abs(ww0 - we) > 15) replace = false;
297  }
298  }
299  if (lclu1 < lclu) {
300  if (c1 != -1 && !deltaraylike && enoughhits) {
301  lclu2 = lclu1;
302  c2 = c1;
303  wb2 = wb1;
304  we2 = we1;
305  tt2 = tt1;
306  dtdw2 = dtdw1;
307  }
308  lclu1 = lclu;
309  c1 = Cls[i][j];
310  wb1 = clusters[Cls[i][j]]->StartWire();
311  we1 = clusters[Cls[i][j]]->EndWire();
312  tt1 = clusters[Cls[i][j]]->StartTick();
313  if (wb1 > we1) {
314  wb1 = clusters[Cls[i][j]]->EndWire();
315  we1 = clusters[Cls[i][j]]->StartWire();
316  tt1 = clusters[Cls[i][j]]->EndTick();
317  }
318  dtdw1 = dtdwstart[Cls[i][j]];
319  }
320  else if (lclu2 < lclu) {
321  if (!deltaraylike && enoughhits && replace) {
322  lclu2 = lclu;
323  c2 = Cls[i][j];
324  wb2 = clusters[Cls[i][j]]->StartWire();
325  we2 = clusters[Cls[i][j]]->EndWire();
326  tt2 = clusters[Cls[i][j]]->StartTick();
327  dtdw2 = dtdwstart[Cls[i][j]];
328  }
329  }
330  }
331  if (c1 != -1 && c2 != -1) {
332  cluvtx[i].push_back(c1);
333  cluvtx[i].push_back(c2);
334 
335  double w1 = clusters[c1]->StartWire();
336  double t1 = clusters[c1]->StartTick();
337  if (clusters[c1]->StartWire() > clusters[c1]->EndWire()) {
338  w1 = clusters[c1]->EndWire();
339  t1 = clusters[c1]->EndTick();
340  }
341  double k1 = dtdwstart[c1];
342  double w2 = clusters[c2]->StartWire();
343  double t2 = clusters[c2]->StartTick();
344  if (clusters[c2]->StartWire() > clusters[c2]->EndWire()) {
345  w1 = clusters[c2]->EndWire();
346  t1 = clusters[c2]->EndTick();
347  }
348  double k2 = dtdwstart[c2];
349  // calculate the vertex
350  if (std::abs(k1 - k2) < 0.5) {
351  vtx_w.push_back(w1);
352  vtx_t.push_back(t1);
353  }
354  else {
355  double t0 = (k1 * k2 * (w1 - w2) + k1 * t2 - k2 * t1) / (k1 - k2);
356  double w0 = (t2 - t1 + k1 * w1 - k2 * w2) / (k1 - k2);
357  vtx_w.push_back(w0);
358  vtx_t.push_back(t0);
359  }
360  }
361  else if (Cls[i].size() >= 1) {
362  if (c1 != -1) {
363  cluvtx[i].push_back(c1);
364  vtx_w.push_back(wb1);
365  vtx_t.push_back(tt1);
366  }
367  else {
368  cluvtx[i].push_back(Cls[i][0]);
369  vtx_w.push_back(clusters[Cls[i][0]]->StartWire());
370  vtx_t.push_back(clusters[Cls[i][0]]->StartTick());
371  }
372  }
373  //save 2D vertex
374  // make an empty art::PtrVector of hits
375  /// \todo should really get the actual vector of hits corresponding to end point
376  /// \todo for now will get all hits from the current cluster
377  std::vector<art::Ptr<recob::Hit>> hits = fmh.at(Cls[i][0]);
378  double totalQ = 0.;
379  for (size_t h = 0; h < hits.size(); ++h)
380  totalQ += hits[h]->Integral();
381 
382  geo::WireID wireID(hits[0]->WireID().Cryostat,
383  hits[0]->WireID().TPC,
384  hits[0]->WireID().Plane,
385  (unsigned int)vtx_w.back()); //for update to EndPoint2D ... WK 4/22/13
386 
387  recob::EndPoint2D vertex(vtx_t.back(),
388  wireID, //for update to EndPoint2D ... WK 4/22/13
389  1,
390  epcol->size(),
391  clusters[Cls[i][0]]->View(),
392  totalQ);
393  epcol->push_back(vertex);
394 
395  util::CreateAssn(evt, *epcol, hits, *assnep);
396  }
397  else {
398  //no cluster found
399  vtx_w.push_back(-1);
400  vtx_t.push_back(-1);
401  }
402  }
403  //std::cout<<vtx_w[0]<<" "<<vtx_t[0]<<" "<<vtx_w[1]<<" "<<vtx_t[1]<<std::endl;
404 
405  Double_t vtxcoord[3];
406  if (Cls[0].size() > 0 && Cls[1].size() > 0) { //ignore w view
407  double Iw0 = (vtx_w[0] + 3.95) * wire_pitch;
408  double Cw0 = (vtx_w[1] + 1.84) * wire_pitch;
409 
410  double It0 = vtx_t[0] - presamplings;
411  It0 *= timepitch;
412  double Ct0 = vtx_t[1] - presamplings;
413  Ct0 *= timepitch;
414  vtxcoord[0] = detProp.ConvertTicksToX(vtx_t[1], 1, 0, 0);
415  vtxcoord[1] = (Cw0 - Iw0) / (2. * std::sin(Angle));
416  vtxcoord[2] = (Cw0 + Iw0) / (2. * std::cos(Angle)) - YC / 2. * std::tan(Angle);
417 
418  double yy, zz;
419  if (vtx_w[0] >= 0 && vtx_w[0] <= 239 && vtx_w[1] >= 0 && vtx_w[1] <= 239) {
420  if (geom->ChannelsIntersect(geom->PlaneWireToChannel(0, (int)((Iw0 / wire_pitch) - 3.95)),
421  geom->PlaneWireToChannel(1, (int)((Cw0 / wire_pitch) - 1.84)),
422  yy,
423  zz)) {
424  //channelsintersect provides a slightly more accurate set of y and z coordinates.
425  // use channelsintersect in case the wires in question do cross.
426  vtxcoord[1] = yy;
427  vtxcoord[2] = zz;
428  }
429  else {
430  vtxcoord[0] = -99999;
431  vtxcoord[1] = -99999;
432  vtxcoord[2] = -99999;
433  }
434  }
435  dtIC->Fill(It0 - Ct0);
436  }
437  else {
438  vtxcoord[0] = -99999;
439  vtxcoord[1] = -99999;
440  vtxcoord[2] = -99999;
441  }
442 
443  /// \todo need to actually make tracks and showers to go into 3D vertex
444  /// \todo currently just passing empty collections to the ctor
445  art::PtrVector<recob::Track> vTracks_vec;
446  art::PtrVector<recob::Shower> vShowers_vec;
447 
448  recob::Vertex the3Dvertex(vtxcoord, vcol->size());
449  vcol->push_back(the3Dvertex);
450 
451  if (vShowers_vec.size() > 0) { util::CreateAssn(evt, *vcol, vShowers_vec, *assnsh); }
452 
453  if (vTracks_vec.size() > 0) { util::CreateAssn(evt, *vcol, vTracks_vec, *assntr); }
454 
455  MF_LOG_VERBATIM("Summary") << std::setfill('-') << std::setw(175) << "-" << std::setfill(' ');
456  MF_LOG_VERBATIM("Summary") << "VertexFinder2D Summary:";
457  for (size_t i = 0; i < epcol->size(); ++i)
458  MF_LOG_VERBATIM("Summary") << epcol->at(i);
459  for (size_t i = 0; i < vcol->size(); ++i)
460  MF_LOG_VERBATIM("Summary") << vcol->at(i);
461 
462  evt.put(std::move(epcol));
463  evt.put(std::move(vcol));
464  evt.put(std::move(assnep));
465  evt.put(std::move(assntr));
466  evt.put(std::move(assnsh));
467  evt.put(std::move(assnh));
468 
469  } // end of produce
470 } // end of vertex namespace
471 
472 // //-----------------------------------------------------------------------------
473 
474 namespace vertex {
475 
477 
478 } // end of vertex namespace
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
code to link reconstructed objects back to the MC truth information
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
AdcChannelData::View View
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
std::string string
Definition: nybbler.cc:12
Planes which measure V.
Definition: geo_types.h:130
std::set< geo::View_t > const & Views() const
Returns a list of possible views in the detector.
constexpr T pow(T x)
Definition: pow.h:72
bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double &y, double &z) const
Returns an intersection point of two channels.
Planes which measure Z direction.
Definition: geo_types.h:132
Definition of vertex object for LArSoft.
Definition: Vertex.h:35
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:250
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
art framework interface to geometry description
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
Definition: type_traits.h:61
T abs(T value)
Planes which measure U.
Definition: geo_types.h:129
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void beginJob()
Definition: Breakpoints.cc:14
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
std::void_t< T > n
def move(depos, offset)
Definition: depos.py:107
double alpha
Definition: doAna.cpp:15
void produce(art::Event &evt) override
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
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.
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
Encapsulate the geometry of a wire.
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Declaration of signal hit object.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
Encapsulate the construction of a single detector plane.
#define MF_LOG_VERBATIM(category)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Provides recob::Track data product.
int trigger_offset(DetectorClocksData const &data)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
TCEvent evt
Definition: DataStructs.cxx:7
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
recob::tracking::Plane Plane
Definition: TrackState.h:17
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:230
Q_EXPORT QTSManip setfill(int f)
Definition: qtextstream.h:337
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.
vertex reconstruction