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 (std::isinf(l.
m) || std::isnan(l.
m) || std::isinf(l.
c) || std::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,
338 pts[0].emplace_back(xpos, r0.z(), 0,
energy);
343 auto const unit = (r1 - r0).Unit();
346 if (perp.z() < 0) perp *= -1;
351 if (dirU.Mag2() == 0) { dirU = perp; }
352 else if (dirV.Mag2() == 0 && fabs(dirU.Dot(perp)) < 0.99) {
358 if (fabs(dirU.Dot(perp)) > 0.99) { pts[1].emplace_back(xpos, r0.Dot(dirU), 1,
energy); }
360 pts[2].emplace_back(xpos, r0.Dot(dirV), 2,
energy);
364 dirs = {dirZ, dirU, dirV};
366 std::default_random_engine gen;
376 const std::vector<recob::Hit>&
hits,
380 if (hits.empty())
return false;
382 std::vector<std::vector<Pt2D>> pts;
383 std::vector<recob::tracking::Vector_t> dirs;
385 GetPts2D(detProp, hits, pts, dirs, wireReadoutGeom);
389 double minz[3] = {+1e9, +1e9, +1e9};
390 double maxz[3] = {-1e9, -1e9, -1e9};
393 minx = std::min(minx, p.x);
394 maxx = std::max(maxx, p.x);
395 minz[
view] = std::min(minz[view], p.z);
396 maxz[
view] = std::max(maxz[view], p.z);
410 std::vector<float> zs;
411 zs.reserve(pts[0].
size());
412 for (
const Pt2D& p : pts[0])
414 auto mid = zs.begin() + zs.size() / 4;
415 if (mid != zs.end()) {
416 std::nth_element(zs.begin(), mid, zs.end());
420 std::vector<HeatMap> hms;
425 std::vector<Line2D> lines;
428 if (lines.empty())
return false;
431 hms.emplace_back(maxz[view] - minz[view], minz[view], maxz[view], maxx - minx, minx, maxx);
437 std::vector<HeatMap> hms_zoom;
440 const double x0 = vtx.X();
441 const double z0 = vtx.Dot(dirs[
view]);
443 std::vector<Line2D> lines;
446 if (lines.empty())
return false;
449 hms_zoom.emplace_back(50, z0 - 2.5, z0 + 2.5, 50, x0 - 2.5, x0 + 2.5);
457 art::TFileDirectory evt_dir =
461 art::TFileDirectory view_dir = evt_dir.mkdir(TString::Format(
"view%d",
view).Data());
463 TGraph* gpts = view_dir.makeAndRegister<TGraph>(
"hits",
"");
465 gpts->SetPoint(gpts->GetN(), p.z, p.x);
467 view_dir.makeAndRegister<TH2F>(
"hmap",
"", *hms[
view].AsTH2());
469 view_dir.makeAndRegister<TH2F>(
"hmap_zoom",
"", *hms_zoom[
view].AsTH2());
471 const double x = vtx.X();
472 const double z = vtx.Dot(dirs[view]);
473 view_dir.makeAndRegister<TGraph>(
"vtx3d",
"", 1, &
z, &
x);
483 auto vtxcol = std::make_unique<std::vector<recob::Vertex>>();
491 if (FindVtx(detProp, *hits, vtx, evt.
event())) {
495 evt.
put(std::move(vtxcol));
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::WireReadoutGeom *wireReadoutGeom)
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.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
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)
Interface for a class providing readout channel mapping to geometry.
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.
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)
View_t View(raw::ChannelID_t const channel) const
Returns the view (wire orientation) on the specified TPC channel.
bool operator<(const Line2D &l) const
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::WireReadoutGeom * wireReadoutGeom
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...
void WireEndPoints(WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
void produce(art::Event &evt) override
Line2D(const Pt2D &a, const Pt2D &b)