129 auto const det_prop =
133 std::vector<art::Ptr<recob::Hit>>
hit;
144 unsigned int numberwires = 0;
145 const unsigned int numbertimesamples = det_prop.ReadOutWindowSize();
146 const float TicksPerBin = numbertimesamples /
fTimeBins;
148 double MatrixAAsum = 0.;
149 double MatrixBBsum = 0.;
150 double MatrixCCsum = 0.;
151 std::vector<double> Cornerness2;
154 double wx[49] = {0.};
155 double wy[49] = {0.};
157 for (
int i = -3; i < 4; ++i) {
158 for (
int j = 3; j > -4; --j) {
166 unsigned int wire = 0;
167 unsigned int wire2 = 0;
168 for (
auto view : geom->
Views()) {
173 while (clusterIter != clusIn.
end()) {
174 if ((*clusterIter)->View() == view) {
175 hit = fmh.at(cinctr);
181 if (hit.size() == 0)
continue;
185 mf::LogInfo(
"EndPointAlg") <<
" --- endpoints check " << numberwires <<
" " << numbertimesamples
188 std::vector<std::vector<double>> MatrixAsum(numberwires);
189 std::vector<std::vector<double>> MatrixBsum(numberwires);
190 std::vector<std::vector<double>> hit_map(numberwires);
193 std::vector<std::vector<int>> hit_loc(numberwires);
195 std::vector<std::vector<double>> Cornerness(numberwires);
197 for (
unsigned int wi = 0; wi < numberwires; ++wi) {
198 hit_map[wi].resize(fTimeBins, 0);
199 hit_loc[wi].resize(fTimeBins, -1);
200 Cornerness[wi].resize(fTimeBins, 0);
201 MatrixAsum[wi].resize(fTimeBins, 0);
202 MatrixBsum[wi].resize(fTimeBins, 0);
204 for (
unsigned int i = 0; i < hit.size(); ++i) {
205 wire = hit[i]->WireID().Wire;
206 const float center = hit[i]->PeakTime(), sigma = hit[i]->RMS();
207 const int iFirstBin = int((center - 3 * sigma) / TicksPerBin),
208 iLastBin = int((center + 3 * sigma) / TicksPerBin);
209 for (
int iBin = iFirstBin; iBin <= iLastBin; ++iBin) {
210 const float bin_center = iBin * TicksPerBin;
211 hit_map[wire][iBin] +=
Gaussian(bin_center, center, sigma);
216 for (
unsigned int wire = 1; wire < numberwires - 1; ++wire) {
218 for (
int timebin = 1; timebin < fTimeBins - 1; ++timebin) {
219 MatrixAsum[wire][timebin] = 0.;
220 MatrixBsum[wire][timebin] = 0.;
222 for (
int i = -3; i <= 3; ++i) {
224 if (windex < 0) windex = 0;
226 else if ((
unsigned int)windex >= numberwires)
227 windex = numberwires - 1;
229 for (
int j = -3; j <= 3; ++j) {
230 tindex = timebin + j;
233 else if (tindex >= fTimeBins)
234 tindex = fTimeBins - 1;
236 MatrixAsum[wire][timebin] += wx[
n] * hit_map[windex][tindex];
237 MatrixBsum[wire][timebin] += wy[
n] * hit_map[windex][tindex];
245 for (
unsigned int wire = 1; wire < numberwires - 1; ++wire) {
247 for (
int timebin = 1; timebin < fTimeBins - 1; ++timebin) {
253 for (
int i = -3; i <= 3; ++i) {
255 if (windex < 0) windex = 0;
257 else if ((
unsigned int)windex >= numberwires)
258 windex = numberwires - 1;
260 for (
int j = -3; j <= 3; ++j) {
261 tindex = timebin + j;
264 else if (tindex >= fTimeBins)
265 tindex = fTimeBins - 1;
267 MatrixAAsum += w[
n] * pow(MatrixAsum[windex][tindex], 2);
268 MatrixBBsum += w[
n] * pow(MatrixBsum[windex][tindex], 2);
269 MatrixCCsum += w[
n] * MatrixAsum[windex][tindex] * MatrixBsum[windex][tindex];
274 if ((MatrixAAsum + MatrixBBsum) > 0)
275 Cornerness[wire][timebin] =
276 (MatrixAAsum * MatrixBBsum - pow(MatrixCCsum, 2)) / (MatrixAAsum + MatrixBBsum);
278 Cornerness[wire][timebin] = 0;
280 if (Cornerness[wire][timebin] > 0) {
281 for (
unsigned int i = 0; i < hit.size(); ++i) {
282 wire2 = hit[i]->WireID().Wire;
284 if (wire == wire2 &&
std::abs(hit[i]->TimeDistanceAsRMS(
285 timebin * (numbertimesamples / fTimeBins))) < 1.) {
287 hit_loc[wire][timebin] = i;
288 Cornerness2.push_back(Cornerness[wire][timebin]);
296 std::sort(Cornerness2.rbegin(), Cornerness2.rend());
298 for (
int vertexnum = 0; vertexnum <
fMaxCorners; ++vertexnum) {
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) {
308 if (Cornerness2.size())
309 if (Cornerness[wire][timebin] < (
fThreshold * Cornerness2[0]))
310 vertexnum = fMaxCorners;
311 vHits.
push_back(hit[hit_loc[wire][timebin]]);
315 for (
size_t vh = 0; vh < vHits.
size(); ++vh)
316 totalQ += vHits[vh]->Integral();
319 hit[hit_loc[wire][timebin]]->
WireID(),
320 Cornerness[wire][timebin],
324 vtxcol.push_back(endpoint);
325 vtxHitsOut.push_back(vHits);
333 double drifttick = det_prop.DriftVelocity(det_prop.Efield(), det_prop.Temperature());
336 double corrfactor = drifttick / wirepitch;
340 (
int)((
fWindow * (numbertimesamples / fTimeBins) * corrfactor) + .5);
342 (int)wire + (
int)((
fWindow * (numbertimesamples /
fTimeBins) * corrfactor) + .5);
344 for (
int timebinout = timebin -
fWindow; timebinout <= timebin +
fWindow;
346 if (std::sqrt(pow(wire - wireout, 2) + pow(timebin - timebinout, 2)) <
348 Cornerness[wireout][timebinout] = 0;
357 if (clusterIter != clusIn.
end()) clusterIter++;
360 unsigned char* outPix =
new unsigned char[fTimeBins * numberwires];
367 for (
unsigned int x = 0;
x < numberwires; ++
x) {
368 cell = (int)(hit_map[
x][
y] * 1000);
369 if (cell > maxCell) { maxCell = cell; }
373 for (
unsigned int x = 0;
x < numberwires; ++
x) {
375 if (maxCell > 0) pix = (int)((1500000 * hit_map[
x][
y]) / maxCell);
376 outPix[y * numberwires +
x] = pix;
382 for (
unsigned int ii = 0; ii < vtxcol.size(); ++ii) {
383 if (vtxcol[ii].View() == (
unsigned int)view) {
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) *
393 vtxcol[ii].
WireID().Wire] = pix;
394 outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
396 vtxcol[ii].
WireID().Wire] = pix;
397 outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
399 vtxcol[ii].
WireID().Wire - 1] = pix;
400 outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
402 vtxcol[ii].
WireID().Wire - 2] = pix;
403 outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
405 vtxcol[ii].
WireID().Wire + 1] = pix;
406 outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
408 vtxcol[ii].
WireID().Wire + 2] = pix;
420 return vtxcol.size();
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
constexpr auto abs(T v)
Returns the absolute value of the argument.
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
void VSSaveBMPFile(const char *fileName, unsigned char *pix, int dx, int dy) const
typename data_t::const_iterator const_iterator
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void push_back(Ptr< U > const &p)
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
PlaneID_t Plane
Index of the plane within its TPC.
double GaussianDerivativeY(int x, int y) const
Detector simulation of raw signals on wires.
double Gaussian(int x, int y, double sigma) const
double GaussianDerivativeX(int x, int y) const
std::set< View_t > const & Views() const
Returns a list of possible views in the detector.
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.