Util.cxx

Go to the documentation of this file.

#include "LEM/Util.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "TVector3.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Cluster.h"

namespace lem
{
  //......................................................................
  int MeanCell(geo::View_t view, const rb::Cluster& slice)
  {
    return (slice.MinCell(view)+slice.MaxCell(view))/2;
  }

  //......................................................................
  void DefaultVertex(const rb::Cluster& slice,
                     int& plane, int& cell, int& cellOtherView,
                     bool reverse)
  {
    art::ServiceHandle<geo::Geometry> geom;

    // Fallback when everything goes wrong.

    // Just put the vertex at the top of the slice in the middle

    plane = reverse ? slice.MaxPlane() : slice.MinPlane();
    if(plane < 0 || plane >= int(geom->NPlanes())) plane = 0; // arghhh
    const geo::View_t view = geom->Plane(plane)->View();
    cell = MeanCell(view, slice);
    cellOtherView = MeanCell(geo::View_t(1-view), slice);
  }

  //......................................................................
  geo::CellUniqueId FindCellRobust(double x, double y, double z)
  {
    art::ServiceHandle<geo::Geometry> geom;

    for(int trial = 0; trial < 4; ++trial){
      // Figure out what cell the vertex was in. Step 3cm in z in case of
      // failure. If that didn't work, it must be in plastic between cells,
      // move laterally in x or y.
      double offsetX = 0, offsetY = 0;
      if(trial == 1) offsetX = -2;
      if(trial == 2) offsetY = -2;
      if(trial == 3) offsetX = offsetY = -2;

      const geo::CellUniqueId id = geom->CellId(x+offsetX, y+offsetY, z,
                                                0, 0, 1, 3);
      // We found one that works
      if(id) return id;
    }
    return 0;
  }

  //......................................................................
  void VertexToPlaneAndCell(const TVector3 vtx, const rb::Cluster& slice,
                            int& plane, int& cell, int& cellOtherView,
                            bool reverse)
  {
    // Picked up by failure cases below
    DefaultVertex(slice, plane, cell, cellOtherView, reverse);

    art::ServiceHandle<geo::Geometry> geom;

    const geo::CellUniqueId id = FindCellRobust(vtx.X(),
                                                vtx.Y(),
                                                vtx.Z());

    // Take defaults
    if(!id) return;

    geom->IdToCell(id, &plane, &cell);

    const int nextPlane = geom->NextPlaneOtherView(plane);
    // Take default in other view

    if(nextPlane == geo::kPLANE_NOT_FOUND) return;

    double xyz[3];
    int ncells = geom->Plane(nextPlane)->Ncells();
    if(cellOtherView >= ncells)
      cellOtherView = ncells-1;

    geom->Plane(nextPlane)->Cell(cellOtherView)->GetCenter(xyz);

    // z position of plane in other view
    double zOther = xyz[2];

    // Sometimes this goes wrong and we have to iterate a bit
    while(true){
      const geo::CellUniqueId idOther = FindCellRobust(vtx.X(),
                                                       vtx.Y(),
                                                       zOther);
      // Take default in other view
      if(!idOther) return;

      // The cell we found must be in the plane we wanted
      int checkPlane;
      geom->IdToCell(idOther, &checkPlane, &cellOtherView);
      // If so, we're good and return the answer
      if(checkPlane == nextPlane) return;
      // Otherwise, walk 1cm at a time in the right direction
      if(checkPlane < nextPlane) zOther += 1;
      if(checkPlane > nextPlane) zOther -= 1;
    }
  }
} // namespace