16 #include "art_root_io/TFileService.h" 18 #include "cetlib/pow.h" 37 Pt2D(
double _x,
double _z,
int _view,
double _energy)
51 : m((b.
x - a.
x) / (b.
z - a.
z))
54 minz(
std::min(a.
z, b.
z))
55 , maxz(
std::max(a.
z, b.
z))
75 const std::vector<recob::Hit>&
hits,
91 , fHitLabel(pset.
get<
std::
string>("HitLabel"))
92 , fSavePlots(pset.
get<
bool>("SavePlots"))
94 produces<std::vector<recob::Vertex>>();
111 const float A = 1 + cet::square(m);
112 const float B = 2 * m * c;
113 const float C = cet::square(c) - cet::square(R);
115 double desc = cet::square(B) - 4 * A * C;
117 if (desc < 0)
return false;
121 z1 = (-B - desc) / (2 * A);
122 z2 = (-B + desc) / (2 * A);
133 std::vector<Line2D>& lines,
138 constexpr
size_t kMaxLines = 10 * 1000 * 1000;
140 const size_t product = (pts.size() * (pts.size() - 1)) / 2;
141 const int stride = product / kMaxLines + 1;
143 lines.reserve(std::min(product, kMaxLines));
145 for (
int offset = 0; offset < stride; ++offset) {
146 for (
unsigned int i = 0; i < pts.size(); ++i) {
147 for (
unsigned int j = i + offset + 1; j < pts.size(); j += stride) {
148 const Line2D l(pts[i], pts[j]);
150 if (isinf(l.
m) || isnan(l.
m) || isinf(l.
c) || isnan(l.
c))
continue;
158 lines.push_back(std::move(l));
159 if (lines.size() == kMaxLines)
goto end;
166 lines.shrink_to_fit();
168 mf::LogInfo() <<
"Made " << lines.size() <<
" lines using stride " << stride
169 <<
" to fit under cap of " << kMaxLines << std::endl;
172 std::sort(lines.begin(), lines.end());
178 const float cosCrit = cos(10 * M_PI / 180.);
179 const float dot = 1 + ma * mb;
180 return cet::square(dot) > (1 + cet::square(ma)) * (1 + cet::square(mb)) * cet::square(cosCrit);
187 constexpr
size_t kMaxPts = 10 * 1000 * 1000;
190 unsigned int jmax = 0;
193 for (
unsigned int i = 0; i + 1 < lines.size(); ++i) {
194 const Line2D a = lines[i];
196 j0 = std::max(j0, i + 1);
197 while (j0 < lines.size() &&
CloseAngles(a.
m, lines[j0].m))
199 jmax = std::max(jmax, j0);
200 while (jmax < lines.size() && !
CloseAngles(a.
m, lines[jmax].m))
206 const size_t product = (lines.size() * (lines.size() - 1)) / 2;
207 const int stride = npts / kMaxPts + 1;
209 mf::LogInfo() <<
"Combining lines to points with stride " << stride << std::endl;
211 mf::LogInfo() << npts <<
" cf " << product <<
" ie " << double(npts) / product << std::endl;
216 for (
unsigned int i = 0; i + 1 < lines.size(); ++i) {
217 const Line2D a = lines[i];
219 j0 = std::max(j0, i + 1);
220 while (j0 < lines.size() &&
CloseAngles(a.
m, lines[j0].m))
222 jmax = std::max(jmax, j0);
223 while (jmax < lines.size() && !
CloseAngles(a.
m, lines[jmax].m))
226 for (
unsigned int j = j0; j < jmax; j += stride) {
227 const Line2D& b = lines[j];
230 const float z = (b.
c - a.
c) / (a.
m - b.
m);
231 const float x = a.
m * z + a.
c;
234 if ((z < a.minz || z > a.
maxz) && (z < b.
minz || z > b.
maxz)) {
235 const int iz = hm.
ZToBin(z);
236 const int ix = hm.
XToBin(x);
237 if (iz >= 0 && iz < hm.Nz && ix >= 0 && ix < hm.
Nx) { hm.
map[iz * hm.
Nx + ix] += stride; }
246 const std::vector<recob::tracking::Vector_t>& dirs) noexcept
248 assert(hs.size() == 3);
249 assert(dirs.size() == 3);
251 const int Nx = hs[0].Nx;
254 M(0, 0) = dirs[0].Y();
255 M(0, 1) = dirs[0].Z();
256 M(1, 0) = dirs[1].Y();
257 M(1, 1) = dirs[1].Z();
260 if (M(0, 0) * M(1, 1) - M(1, 0) * M(0, 1) == 0)
return {};
264 float bestscore = -1;
268 std::vector<float> colMax[3];
271 for (
int iz = 0; iz < hs[
view].Nz; ++iz) {
272 colMax[
view][iz] = *std::max_element(&hs[
view].map[Nx * iz], &hs[
view].map[Nx * (iz + 1)]);
276 for (
int iz = 0; iz < hs[0].Nz; ++iz) {
277 const float z = hs[0].ZBinCenter(iz);
278 const float bonus = 1;
280 for (
int iu = 0; iu < hs[1].Nz; ++iu) {
281 const float u = hs[1].ZBinCenter(iu);
286 const TVectorD
r = M * p;
287 const float v = r[0] * dirs[2].Y() + r[1] * dirs[2].Z();
288 const int iv = hs[2].ZToBin(v);
289 if (iv < 0 || iv >= hs[2].Nz)
continue;
290 const double y =
r(0);
293 if (colMax[0][iz] + colMax[1][iu] + colMax[2][iv] < bestscore)
continue;
295 const float* h0 = &hs[0].map[Nx * iz];
296 const float*
h1 = &hs[1].map[Nx * iu];
297 const float*
h2 = &hs[2].map[Nx * iv];
300 for (
int ix = 1; ix < Nx - 1; ++ix) {
301 const float score = bonus * (h0[ix] + h1[ix] + h2[ix]);
303 if (score > bestscore) {
309 if (bestix != -1) { bestr.SetXYZ(hs[0].XBinCenter(bestix), y, z); }
318 const std::vector<recob::Hit>&
hits,
320 std::vector<recob::tracking::Vector_t>& dirs,
339 pts[0].emplace_back(xpos, r0.z(), 0,
energy);
344 auto const unit = (r1 - r0).Unit();
347 if (perp.z() < 0) perp *= -1;
352 if (dirU.Mag2() == 0) { dirU = perp; }
353 else if (dirV.Mag2() == 0 && fabs(dirU.Dot(perp)) < 0.99) {
359 if (fabs(dirU.Dot(perp)) > 0.99) { pts[1].emplace_back(xpos, r0.Dot(dirU), 1,
energy); }
361 pts[2].emplace_back(xpos, r0.Dot(dirV), 2,
energy);
365 dirs = {dirZ, dirU, dirV};
367 std::default_random_engine gen;
377 const std::vector<recob::Hit>&
hits,
381 if (hits.empty())
return false;
383 std::vector<std::vector<Pt2D>> pts;
384 std::vector<recob::tracking::Vector_t> dirs;
386 GetPts2D(detProp, hits, pts, dirs, geom);
390 double minz[3] = {+1e9, +1e9, +1e9};
391 double maxz[3] = {-1e9, -1e9, -1e9};
394 minx = std::min(minx, p.x);
395 maxx = std::max(maxx, p.x);
396 minz[
view] = std::min(minz[view], p.z);
397 maxz[
view] = std::max(maxz[view], p.z);
411 std::vector<float> zs;
412 zs.reserve(pts[0].
size());
413 for (
const Pt2D& p : pts[0])
415 auto mid = zs.begin() + zs.size() / 4;
416 if (mid != zs.end()) {
417 std::nth_element(zs.begin(), mid, zs.end());
421 std::vector<HeatMap> hms;
426 std::vector<Line2D> lines;
429 if (lines.empty())
return false;
432 hms.emplace_back(maxz[view] - minz[view], minz[view], maxz[view], maxx - minx, minx, maxx);
438 std::vector<HeatMap> hms_zoom;
441 const double x0 = vtx.X();
442 const double z0 = vtx.Dot(dirs[
view]);
444 std::vector<Line2D> lines;
447 if (lines.empty())
return false;
450 hms_zoom.emplace_back(50, z0 - 2.5, z0 + 2.5, 50, x0 - 2.5, x0 + 2.5);
458 art::TFileDirectory evt_dir =
462 art::TFileDirectory view_dir = evt_dir.mkdir(TString::Format(
"view%d",
view).Data());
464 TGraph* gpts = view_dir.makeAndRegister<TGraph>(
"hits",
"");
466 gpts->SetPoint(gpts->GetN(), p.z, p.x);
468 view_dir.makeAndRegister<TH2F>(
"hmap",
"", *hms[
view].AsTH2());
470 view_dir.makeAndRegister<TH2F>(
"hmap_zoom",
"", *hms_zoom[
view].AsTH2());
472 const double x = vtx.X();
473 const double z = vtx.Dot(dirs[view]);
474 view_dir.makeAndRegister<TGraph>(
"vtx3d",
"", 1, &
z, &
x);
484 auto vtxcol = std::make_unique<std::vector<recob::Vertex>>();
492 if (FindVtx(detProp, *hits, vtx, evt.
event())) {
496 evt.
put(std::move(vtxcol));
int XToBin(double x) const
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Declaration of signal hit object.
bool FindVtx(const detinfo::DetectorPropertiesData &detProp, const std::vector< recob::Hit > &hits, recob::tracking::Point_t &vtx, int evt) const
tracking::SMatrixSym33 SMatrixSym33
int ZToBin(double z) const
Planes which measure Z direction.
void GetPts2D(const detinfo::DetectorPropertiesData &detProp, const std::vector< recob::Hit > &hits, std::vector< std::vector< Pt2D >> &pts, std::vector< recob::tracking::Vector_t > &dirs, const geo::GeometryCore *geom)
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
tracking::Point_t Point_t
Pt2D(double _x, double _z, int _view, double _energy)
#define DEFINE_ART_MODULE(klass)
void MapFromLines(const std::vector< Line2D > &lines, HeatMap &hm)
bool IntersectsCircle(float m, float c, float z0, float x0, float R, float &z1, float &z2)
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 ...
bool operator<(const Pt2D &p) const
recob::tracking::Point_t FindPeak3D(const std::vector< HeatMap > &hs, const std::vector< recob::tracking::Vector_t > &dirs) noexcept
EventNumber_t event() const
constexpr auto dot(Vector const &a, OtherVector const &b)
Return cross product of two vectors.
Description of geometry of one entire detector.
Detector simulation of raw signals on wires.
double ConvertTicksToX(double ticks, int p, int t, int c) const
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
bool CloseAngles(float ma, float mb)
bool operator<(const Line2D &l) const
View_t View(PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
void LinesFromPoints(const std::vector< Pt2D > &pts, std::vector< Line2D > &lines, float z0=0, float x0=0, float R=-1)
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
2D representation of charge deposited in the TDC/wire plane
const geo::GeometryCore * geom
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...
art framework interface to geometry description
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
void WireEndPoints(WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
void produce(art::Event &evt) override
Line2D(const Pt2D &a, const Pt2D &b)