LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
EndPointAlg.cxx
Go to the documentation of this file.
1 //
3 // EndPointAlg class
4 //
5 // joshua.spitz@yale.edu
6 //
7 // This algorithm is designed to find (weak) vertices from hits after deconvolution and hit finding.
8 // A weak end point is a end point that has been found using a dedicated end point finding algorithm only. A
9 // strong end point is a end point that has been found using a dedicated end point finding algorithm and matched
10 // to a crossing of two or more HoughLineFinder lines. The End PointMatch module finds strong vertices.
17 //Thanks to B. Morgan of U. of Warwick for comments and suggestions
18 
19 #include <algorithm>
20 #include <fstream>
21 #include <math.h>
22 #include <vector>
23 
24 #include "TMath.h"
25 #include "TString.h"
26 
27 // Framework includes
33 #include "cetlib/pow.h"
35 
46 
47 //-----------------------------------------------------------------------------
49 {
50  fTimeBins = p.get<int>("TimeBins");
51  fMaxCorners = p.get<int>("MaxCorners");
52  fGsigma = p.get<double>("Gsigma");
53  fWindow = p.get<int>("Window");
54  fThreshold = p.get<double>("Threshold");
55  fSaveVertexMap = p.get<int>("SaveVertexMap");
56 }
57 
58 //-----------------------------------------------------------------------------
59 double cluster::EndPointAlg::Gaussian(int x, int y, double sigma) const
60 {
61  double const Norm = 1. / std::sqrt(2 * TMath::Pi() * cet::square(sigma));
62  return Norm * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(sigma)));
63 }
64 
65 //-----------------------------------------------------------------------------
67 {
68  double const Norm = 1. / (std::sqrt(2 * TMath::Pi()) * cet::cube(fGsigma));
69  return Norm * (-x) * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(fGsigma)));
70 }
71 
72 //-----------------------------------------------------------------------------
74 {
75  double const Norm = 1. / (std::sqrt(2 * TMath::Pi()) * cet::cube(fGsigma));
76  return Norm * (-y) * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(fGsigma)));
77 }
78 
79 //-----------------------------------------------------------------------------
80 //this method saves a BMP image of the vertex map space, which can be viewed with gimp
81 void cluster::EndPointAlg::VSSaveBMPFile(const char* fileName,
82  unsigned char* pix,
83  int dx,
84  int dy) const
85 {
86  std::ofstream bmpFile(fileName, std::ios::binary);
87  bmpFile.write("B", 1);
88  bmpFile.write("M", 1);
89  int bitsOffset = 54 + 256 * 4;
90  int size = bitsOffset + dx * dy; //header plus 256 entry LUT plus pixels
91  bmpFile.write((const char*)&size, 4);
92  int reserved = 0;
93  bmpFile.write((const char*)&reserved, 4);
94  bmpFile.write((const char*)&bitsOffset, 4);
95  int bmiSize = 40;
96  bmpFile.write((const char*)&bmiSize, 4);
97  bmpFile.write((const char*)&dx, 4);
98  bmpFile.write((const char*)&dy, 4);
99  short planes = 1;
100  bmpFile.write((const char*)&planes, 2);
101  short bitCount = 8;
102  bmpFile.write((const char*)&bitCount, 2);
103  int i, temp = 0;
104  for (i = 0; i < 6; i++)
105  bmpFile.write((const char*)&temp, 4); // zero out optional color info
106  // write a linear LUT
107  char lutEntry[4]; // blue,green,red
108  lutEntry[3] = 0; // reserved part
109  for (i = 0; i < 256; i++) {
110  lutEntry[0] = i;
111  lutEntry[1] = i + 1;
112  lutEntry[2] = i + 2;
113  bmpFile.write(lutEntry, sizeof lutEntry);
114  }
115  // write the actual pixels
116  bmpFile.write((const char*)pix, dx * dy);
117 }
118 
119 //......................................................
121  std::vector<recob::EndPoint2D>& vtxcol,
123  art::Event const& evt,
124  std::string const& label) const
125 {
126  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
127  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
128  auto const det_prop =
130 
131  //Point to a collection of vertices to output.
132  std::vector<art::Ptr<recob::Hit>> hit;
133 
134  art::FindManyP<recob::Hit> fmh(clusIn, evt, label);
135 
136  int flag = 0;
137  int windex = 0; // the wire index to make sure the end point finder does not
138  // fall off the edge of the hit map
139  int tindex = 0; // the time index to make sure the end point finder does not
140  // fall off the edge of the hit map
141  int n = 0; // index of window cell. There are 49 cells in the 7X7 Gaussian and
142  // Gaussian derivative windows
143  unsigned int numberwires = 0;
144  const unsigned int numbertimesamples = det_prop.ReadOutWindowSize();
145  const float TicksPerBin = numbertimesamples / fTimeBins;
146 
147  double MatrixAAsum = 0.;
148  double MatrixBBsum = 0.;
149  double MatrixCCsum = 0.;
150  std::vector<double> Cornerness2;
151  //gaussian window definitions. The cell weights are calculated here to help the algorithm's speed
152  double w[49] = {0.};
153  double wx[49] = {0.};
154  double wy[49] = {0.};
155  int ctr = 0;
156  for (int i = -3; i < 4; ++i) {
157  for (int j = 3; j > -4; --j) {
158  w[ctr] = Gaussian(i, j, fGsigma);
159  wx[ctr] = GaussianDerivativeX(i, j);
160  wy[ctr] = GaussianDerivativeY(i, j);
161  ++ctr;
162  }
163  }
164 
165  unsigned int wire = 0;
166  unsigned int wire2 = 0;
167  for (auto view : wireReadoutGeom.Views()) {
170  hit.clear();
171  size_t cinctr = 0;
172  while (clusterIter != clusIn.end()) {
173  if ((*clusterIter)->View() == view) {
174  hit = fmh.at(cinctr);
175  } //end if cluster is in the correct view
176  clusterIter++;
177  ++cinctr;
178  }
179 
180  if (hit.size() == 0) continue;
181 
182  geo::WireID wid = hit[0]->WireID();
183  numberwires = wireReadoutGeom.Plane(wid.asPlaneID()).Nwires();
184  mf::LogInfo("EndPointAlg") << " --- endpoints check " << numberwires << " " << numbertimesamples
185  << " " << fTimeBins;
186 
187  std::vector<std::vector<double>> MatrixAsum(numberwires);
188  std::vector<std::vector<double>> MatrixBsum(numberwires);
189  std::vector<std::vector<double>> hit_map(numberwires); //the map of hits
190 
191  //the index of the hit that corresponds to the potential corner
192  std::vector<std::vector<int>> hit_loc(numberwires);
193 
194  std::vector<std::vector<double>> Cornerness(numberwires); //the "weight" of a corner
195 
196  for (unsigned int wi = 0; wi < numberwires; ++wi) {
197  hit_map[wi].resize(fTimeBins, 0);
198  hit_loc[wi].resize(fTimeBins, -1);
199  Cornerness[wi].resize(fTimeBins, 0);
200  MatrixAsum[wi].resize(fTimeBins, 0);
201  MatrixBsum[wi].resize(fTimeBins, 0);
202  }
203  for (unsigned int i = 0; i < hit.size(); ++i) {
204  wire = hit[i]->WireID().Wire;
205  const float center = hit[i]->PeakTime(), sigma = hit[i]->RMS();
206  const int iFirstBin = int((center - 3 * sigma) / TicksPerBin),
207  iLastBin = int((center + 3 * sigma) / TicksPerBin);
208  for (int iBin = iFirstBin; iBin <= iLastBin; ++iBin) {
209  const float bin_center = iBin * TicksPerBin;
210  hit_map[wire][iBin] += Gaussian(bin_center, center, sigma);
211  }
212  }
213 
214  // Gaussian derivative convolution
215  for (unsigned int wire = 1; wire < numberwires - 1; ++wire) {
216 
217  for (int timebin = 1; timebin < fTimeBins - 1; ++timebin) {
218  MatrixAsum[wire][timebin] = 0.;
219  MatrixBsum[wire][timebin] = 0.;
220  n = 0;
221  for (int i = -3; i <= 3; ++i) {
222  windex = wire + i;
223  if (windex < 0) windex = 0;
224  // this is ok, because the line before makes sure it's not negative
225  else if ((unsigned int)windex >= numberwires)
226  windex = numberwires - 1;
227 
228  for (int j = -3; j <= 3; ++j) {
229  tindex = timebin + j;
230  if (tindex < 0)
231  tindex = 0;
232  else if (tindex >= fTimeBins)
233  tindex = fTimeBins - 1;
234 
235  MatrixAsum[wire][timebin] += wx[n] * hit_map[windex][tindex];
236  MatrixBsum[wire][timebin] += wy[n] * hit_map[windex][tindex];
237  ++n;
238  } // end loop over j
239  } // end loop over i
240  } // end loop over time bins
241  } // end loop over wires
242 
243  //calculate the cornerness of each pixel while making sure not to fall off the hit map.
244  for (unsigned int wire = 1; wire < numberwires - 1; ++wire) {
245 
246  for (int timebin = 1; timebin < fTimeBins - 1; ++timebin) {
247  MatrixAAsum = 0.;
248  MatrixBBsum = 0.;
249  MatrixCCsum = 0.;
250  //Gaussian smoothing convolution
251  n = 0;
252  for (int i = -3; i <= 3; ++i) {
253  windex = wire + i;
254  if (windex < 0) windex = 0;
255  // this is ok, because the line before makes sure it's not negative
256  else if ((unsigned int)windex >= numberwires)
257  windex = numberwires - 1;
258 
259  for (int j = -3; j <= 3; ++j) {
260  tindex = timebin + j;
261  if (tindex < 0)
262  tindex = 0;
263  else if (tindex >= fTimeBins)
264  tindex = fTimeBins - 1;
265 
266  MatrixAAsum += w[n] * pow(MatrixAsum[windex][tindex], 2);
267  MatrixBBsum += w[n] * pow(MatrixBsum[windex][tindex], 2);
268  MatrixCCsum += w[n] * MatrixAsum[windex][tindex] * MatrixBsum[windex][tindex];
269  ++n;
270  } // end loop over j
271  } // end loop over i
272 
273  if ((MatrixAAsum + MatrixBBsum) > 0)
274  Cornerness[wire][timebin] =
275  (MatrixAAsum * MatrixBBsum - pow(MatrixCCsum, 2)) / (MatrixAAsum + MatrixBBsum);
276  else
277  Cornerness[wire][timebin] = 0;
278 
279  if (Cornerness[wire][timebin] > 0) {
280  for (unsigned int i = 0; i < hit.size(); ++i) {
281  wire2 = hit[i]->WireID().Wire;
282  //make sure the end point candidate coincides with an actual hit.
283  if (wire == wire2 && std::abs(hit[i]->TimeDistanceAsRMS(
284  timebin * (numbertimesamples / fTimeBins))) < 1.) {
285  //this index keeps track of the hit number
286  hit_loc[wire][timebin] = i;
287  Cornerness2.push_back(Cornerness[wire][timebin]);
288  break;
289  }
290  } // end loop over hits
291  } // end if cornerness > 0
292  } // end loop over time bins
293  } // end wire loop
294 
295  std::sort(Cornerness2.rbegin(), Cornerness2.rend());
296 
297  auto const& plane = wireReadoutGeom.Plane({0, 0, 0});
298  for (int vertexnum = 0; vertexnum < fMaxCorners; ++vertexnum) {
299  flag = 0;
300  for (unsigned int wire = 0; wire < numberwires && flag == 0; ++wire) {
301  for (int timebin = 0; timebin < fTimeBins && flag == 0; ++timebin) {
302  if (Cornerness2.size() > (unsigned int)vertexnum)
303  if (Cornerness[wire][timebin] == Cornerness2[vertexnum] &&
304  Cornerness[wire][timebin] > 0. && hit_loc[wire][timebin] > -1) {
305  ++flag;
306 
307  //thresholding
308  if (Cornerness2.size())
309  if (Cornerness[wire][timebin] < (fThreshold * Cornerness2[0]))
310  vertexnum = fMaxCorners;
311  vHits.push_back(hit[hit_loc[wire][timebin]]);
312 
313  // get the total charge from the associated hits
314  double totalQ = 0.;
315  for (size_t vh = 0; vh < vHits.size(); ++vh)
316  totalQ += vHits[vh]->Integral();
317 
318  recob::EndPoint2D endpoint(hit[hit_loc[wire][timebin]]->PeakTime(),
319  hit[hit_loc[wire][timebin]]->WireID(),
320  Cornerness[wire][timebin],
321  vtxcol.size(),
322  view,
323  totalQ);
324  vtxcol.push_back(endpoint);
325  vtxHitsOut.push_back(vHits);
326  vHits.clear();
327 
328  // non-maximal suppression on a square window. The wire coordinate units are
329  // converted to time ticks so that the window is truly square.
330  // Note that there are 1/0.0743=13.46 time samples per 4.0 mm (wire pitch in ArgoNeuT),
331  // assuming a 1.5 mm/us drift velocity for a 500 V/cm E-field
332 
333  double drifttick = det_prop.DriftVelocity(det_prop.Efield(), det_prop.Temperature());
334  drifttick *= sampling_rate(clock_data) * 1.e-3;
335  double wirepitch = plane.WirePitch();
336  double corrfactor = drifttick / wirepitch;
337 
338  for (int wireout =
339  (int)wire -
340  (int)((fWindow * (numbertimesamples / fTimeBins) * corrfactor) + .5);
341  wireout <=
342  (int)wire + (int)((fWindow * (numbertimesamples / fTimeBins) * corrfactor) + .5);
343  ++wireout) {
344  for (int timebinout = timebin - fWindow; timebinout <= timebin + fWindow;
345  timebinout++) {
346  if (std::sqrt(pow(wire - wireout, 2) + pow(timebin - timebinout, 2)) <
347  fWindow) //circular window
348  Cornerness[wireout][timebinout] = 0;
349  }
350  }
351  }
352  }
353  }
354  }
355  Cornerness2.clear();
356  hit.clear();
357  if (clusterIter != clusIn.end()) clusterIter++;
358 
359  if ((int)wid.Plane == fSaveVertexMap) {
360  unsigned char* outPix = new unsigned char[fTimeBins * numberwires];
361  //finds the maximum cell in the map for image scaling
362  int cell = 0;
363  int pix = 0;
364  int maxCell = 0;
365  //int xmaxx, ymaxx;
366  for (int y = 0; y < fTimeBins; ++y) {
367  for (unsigned int x = 0; x < numberwires; ++x) {
368  cell = (int)(hit_map[x][y] * 1000);
369  if (cell > maxCell) { maxCell = cell; }
370  }
371  }
372  for (int y = 0; y < fTimeBins; ++y) {
373  for (unsigned int x = 0; x < numberwires; ++x) {
374  //scales the pixel weights based on the maximum cell value
375  if (maxCell > 0) pix = (int)((1500000 * hit_map[x][y]) / maxCell);
376  outPix[y * numberwires + x] = pix;
377  }
378  }
379 
380  // add 3x3 pixel squares to with the harris vertex finders to the .bmp file
381 
382  for (unsigned int ii = 0; ii < vtxcol.size(); ++ii) {
383  if (vtxcol[ii].View() == (unsigned int)view) {
384  pix = (int)(255);
385  outPix[(int)(vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) * numberwires +
386  vtxcol[ii].WireID().Wire] = pix;
387  outPix[(int)(vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) * numberwires +
388  vtxcol[ii].WireID().Wire - 1] = pix;
389  outPix[(int)(vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) * numberwires +
390  vtxcol[ii].WireID().Wire + 1] = pix;
391  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
392  numberwires +
393  vtxcol[ii].WireID().Wire] = pix;
394  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
395  numberwires +
396  vtxcol[ii].WireID().Wire] = pix;
397  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
398  numberwires +
399  vtxcol[ii].WireID().Wire - 1] = pix;
400  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
401  numberwires +
402  vtxcol[ii].WireID().Wire - 2] = pix;
403  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
404  numberwires +
405  vtxcol[ii].WireID().Wire + 1] = pix;
406  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
407  numberwires +
408  vtxcol[ii].WireID().Wire + 2] = pix;
409  }
410  }
411 
412  VSSaveBMPFile(Form("harrisvertexmap_%d_%d.bmp", (*clusIn.begin())->ID(), wid.Plane),
413  outPix,
414  numberwires,
415  fTimeBins);
416  delete[] outPix;
417  }
418  } // end loop over views
419 
420  return vtxcol.size();
421 }
Float_t x
Definition: compare.C:6
algorithm to find 2D endpoints
Encapsulate the construction of a single cyostat .
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
iterator begin()
Definition: PtrVector.h:217
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
constexpr auto abs(T v)
Returns the absolute value of the argument.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
size_t EndPoint(const art::PtrVector< recob::Cluster > &clusIn, std::vector< recob::EndPoint2D > &vtxcol, std::vector< art::PtrVector< recob::Hit >> &vtxHitsOut, art::Event const &evt, std::string const &label) const
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
void VSSaveBMPFile(const char *fileName, unsigned char *pix, int dx, int dy) const
Definition: EndPointAlg.cxx:81
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
EndPointAlg(fhicl::ParameterSet const &pset)
Definition: EndPointAlg.cxx:48
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
T get(std::string const &key) const
Definition: ParameterSet.h:314
iterator end()
Definition: PtrVector.h:231
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
Declaration of cluster object.
size_type size() const
Definition: PtrVector.h:302
double GaussianDerivativeY(int x, int y) const
Definition: EndPointAlg.cxx:73
Detector simulation of raw signals on wires.
double Gaussian(int x, int y, double sigma) const
Definition: EndPointAlg.cxx:59
Encapsulate the construction of a single detector plane .
double GaussianDerivativeX(int x, int y) const
Definition: EndPointAlg.cxx:66
Char_t n[5]
TCEvent evt
Definition: DataStructs.cxx:8
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
void clear()
Definition: PtrVector.h:533
Float_t w
Definition: plot.C:20
Encapsulate the construction of a single detector plane .
constexpr PlaneID const & asPlaneID() const
Conversion to WireID (for convenience of notation).
Definition: geo_types.h:466