APAGeometryAlg.cxx

Go to the documentation of this file.

//
// \fileAPAGeometryAlg.cxx
//
// tylerdalion@gmail.com
//
//  Geometry interface to smooth over the awkwardness of fitting an
//  APA into LArSoft. It is geared towards making reconstruction simpler
//

//Framework includes:
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "larcorealg/CoreUtils/NumericUtils.h"
#include "larcorealg/Geometry/CryostatGeo.h"
#include "larcorealg/Geometry/TPCGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larreco/RecoAlg/APAGeometryAlg.h"

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <iostream>

namespace apa {

  APAGeometryAlg::APAGeometryAlg(fhicl::ParameterSet const& pset)
  {
    this->reconfigure(pset);
    this->Init();
  }

  //----------------------------------------------------------
  APAGeometryAlg::APAGeometryAlg()
  {
    this->Init();
  }

  //----------------------------------------------------------
  void APAGeometryAlg::reconfigure(fhicl::ParameterSet const& /*p*/) {}

  //----------------------------------------------------------
  void APAGeometryAlg::Init()
  {

    // find the number of channels per APA
    uint32_t channel = 0;

    // old logic -- segfaults if there is just one APA
    //while( fGeom->ChannelToWire(channel+1)[0].TPC < 2 ) channel++;
    //fChannelsPerAPA = channel + 1;

    fChannelsPerAPA = fGeom->Nchannels();
    for (channel = 0; channel < fGeom->Nchannels(); ++channel) {
      if (fGeom->ChannelToWire(channel)[0].TPC > 1) {
        fChannelsPerAPA = channel;
        break;
      }
    }

    // Step through channel c and find the view boundaries, until
    // outside of first APA - these help optimize ChannelToAPAView
    // (very dependent on the conventions implimented in the channel map)
    fFirstU = 0;
    uint32_t c = 1;
    geo::WireID wid = (fGeom->ChannelToWire(c))[0];
    geo::WireID lastwid;
    while (wid.TPC < 2) {

      if (fGeom->View(c) == geo::kV && fGeom->View(c - 1) == geo::kU) {
        fLastU = c - 1;
        fFirstV = c;
      }

      if (fGeom->View(c) == geo::kZ && fGeom->View(c - 1) == geo::kV) {
        fLastV = c - 1;
        fFirstZ0 = c;
      }

      if (wid.TPC == lastwid.TPC + 1) {
        fLastZ0 = c - 1;
        fFirstZ1 = c;
      }

      lastwid = wid;
      c++;
      if (c >= fGeom->Nchannels()) break;
      wid = (fGeom->ChannelToWire(c))[0]; // for the while condition
    }

    fLastZ1 = c - 1;

    if (fLastZ1 + 1 != fChannelsPerAPA)
      throw cet::exception("APAGeometryAlg") << "Channel boundaries are inconsistent.\n";

    // some other things that will be needed
    fAPAsPerCryo = fGeom->NTPC() / 2;
    fChannelRange[0] = (fLastU - fFirstU + 1) * fGeom->WirePitch(geo::kU);
    fChannelRange[1] = (fLastV - fFirstV + 1) * fGeom->WirePitch(geo::kV);
  }

  //----------------------------------------------------------
  void APAGeometryAlg::ChannelToAPA(uint32_t chan, unsigned int& apa, unsigned int& cryo) const
  {

    cryo = chan / (fAPAsPerCryo * fChannelsPerAPA);

    // Number apa uniquely across cryostats so that
    // apa to recob::Object maps are easy to work with.
    // If we decide to reset APA number per cryo, uncomment:
    //chan -= cryo*fAPAsPerCryo*fChannelsPerAPA;
    apa = chan / fChannelsPerAPA;

    return;
  }

  //----------------------------------------------------------
  unsigned int APAGeometryAlg::ChannelToAPA(uint32_t chan) const
  {
    return chan / fChannelsPerAPA;
  }

  //----------------------------------------------------------
  unsigned int APAGeometryAlg::ChannelsInView(geo::View_t geoview) const
  {

    switch (geoview) {
    default: return 0;
    case geo::kU: return ChannelsInAPAView(kU);
    case geo::kV: return ChannelsInAPAView(kV);
    case geo::kZ:
      if (ChannelsInAPAView(kZ0) == ChannelsInAPAView(kZ1))
        return ChannelsInAPAView(kZ0);
      else
        throw cet::exception("ChannelsInView")
          << "Both Z sides should have the same amount of channels\n";
    }
  }

  //----------------------------------------------------------
  unsigned int APAGeometryAlg::ChannelsInAPAView(APAView_t apaview) const
  {

    switch (apaview) {
    default: return 0;
    case kU: return fFirstV - fFirstU;
    case kV: return fFirstZ0 - fFirstV;
    case kZ0: return fFirstZ1 - fFirstZ0;
    case kZ1: return fLastZ1 - fFirstZ1 + 1;
    }
  }

  //----------------------------------------------------------
  uint32_t APAGeometryAlg::FirstChannelInView(geo::View_t geoview,
                                              unsigned int apa,
                                              unsigned int cryo) const
  {

    switch (geoview) {
    default: return 0 + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
    case geo::kU: return fFirstU + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
    case geo::kV: return fFirstV + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
    case geo::kZ:
      //TODO: would need tpc number for the rest of this
      return fFirstZ0 + (uint32_t)(apa + cryo * fAPAsPerCryo) * fChannelsPerAPA;
    }
  }

  //----------------------------------------------------------
  uint32_t APAGeometryAlg::FirstChannelInView(uint32_t chan) const
  {

    geo::View_t geoview = fGeom->View(chan);
    unsigned int apa, cryo;
    this->ChannelToAPA(chan, apa, cryo);
    return this->FirstChannelInView(geoview, chan);
  }

  //----------------------------------------------------------
  uint32_t APAGeometryAlg::FirstChannelInView(geo::View_t geoview, uint32_t chan) const
  {

    unsigned int apa, cryo;
    this->ChannelToAPA(chan, apa, cryo);
    return this->FirstChannelInView(geoview, apa, cryo);
  }

  //----------------------------------------------------------
  APAView_t APAGeometryAlg::APAView(uint32_t chan) const
  {

    // it seems trivial to do this for U and V, but this gives a side to
    // geo::kZ, unlike Geometry::View(c), as is often needed in disambiguation

    geo::View_t view = fGeom->View(chan);
    switch (view) {
    default: break;
    case geo::kU: return kU;
    case geo::kV: return kV;
    case geo::kZ:
      unsigned int modchan = chan % fChannelsPerAPA;
      // Channel mapping number in the order of U, V, Z0, then Z1
      if (modchan > fLastV && modchan < fFirstZ1) return kZ0;
      if (modchan > fLastZ0) return kZ1;
    }

    return kUnknown;
  }

  //----------------------------------------------------------
  std::vector<geo::WireID> APAGeometryAlg::ChanSegsPerSide(uint32_t chan, unsigned int side) const
  {

    std::vector<geo::WireID> wids = fGeom->ChannelToWire(chan);
    return this->ChanSegsPerSide(wids, side);
  }

  //----------------------------------------------------------
  std::vector<geo::WireID> APAGeometryAlg::ChanSegsPerSide(std::vector<geo::WireID> wids,
                                                           unsigned int side) const
  {
    // Given a vector of wireIDs and an APA side, return
    // the wireIDs the the tpc side where tpc%2 = side

    std::vector<geo::WireID> thisSide;

    for (size_t i = 0; i < wids.size(); i++)
      if (wids[i].TPC % 2 == side) thisSide.push_back(wids[i]);

    return thisSide;
  }

  //----------------------------------------------------------
  geo::WireID APAGeometryAlg::NearestWireIDOnChan(const double WorldLoc[3],
                                                  uint32_t chan,
                                                  unsigned int const plane,
                                                  unsigned int const tpc,
                                                  unsigned int const cstat) const
  {

    std::vector<geo::WireID> cWids = fGeom->ChannelToWire(chan);

    if (cWids[0].Cryostat != cstat)
      throw cet::exception("APAGeometryAlg")
        << "Channel " << chan << "not in cryostat " << cstat << "\n";
    if (std::floor(cWids[0].TPC / 2) != std::floor(tpc / 2))
      throw cet::exception("APAGeometryAlg")
        << "Channel " << chan << "not in APA " << std::floor(tpc / 2) << "\n";

    // special case for vertical wires
    if (fGeom->View(chan) == geo::kZ) return fGeom->ChannelToWire(chan)[0];

    unsigned int xyzWire = fGeom
                             ->NearestWireID(geo::Point_t{WorldLoc[0], WorldLoc[1], WorldLoc[2]},
                                             geo::PlaneID{cstat, tpc, plane})
                             .Wire;

    // The desired wire ID will be the only channel
    // segment within half the channel range.
    geo::WireID wid;
    for (size_t i = 0; i < cWids.size(); i++) {
      if (cWids[i].TPC != tpc) continue;
      if (std::abs((int)cWids[i].Wire - (int)xyzWire) < fChannelRange[plane] / 2) wid = cWids[i];
    }

    return wid;
  }

  //----------------------------------------------------------
  bool APAGeometryAlg::LineSegChanIntersect(TVector3 xyzStart,
                                            TVector3 xyzEnd,
                                            uint32_t chan,
                                            std::vector<geo::WireID>& widsCrossed,
                                            bool ExtendLine = true) const
  {
    // This assumes a smooth wire numbering, and that the line seg is contained in a tpc.
    // Meant for use with the approximate line calculated
    // by matching cluster endpoints in disambiguation.

    // Find tpc, use midpoint in case start/end is on a boundary
    geo::Point_t const xyzMid{
      (xyzStart[0] + xyzEnd[0]) / 2, (xyzStart[1] + xyzEnd[1]) / 2, (xyzStart[2] + xyzEnd[2]) / 2};
    auto const& tpcid = fGeom->PositionToTPCID(xyzMid);

    // Find the nearest wire number to the line segment endpoints
    std::vector<geo::WireID> wids = fGeom->ChannelToWire(chan);
    geo::PlaneID const planeID{tpcid, wids[0].Plane};
    using geo::vect::toPoint;
    unsigned int startW = fGeom->NearestWireID(toPoint(xyzStart), planeID).Wire;
    unsigned int endW = fGeom->NearestWireID(toPoint(xyzEnd), planeID).Wire;

    if (startW > endW) std::swap(startW, endW);

    // Loop through wireIDs and check for intersection, if in the right TPC
    for (size_t w = 0; w < wids.size(); w++) {
      if (wids[w].TPC != tpcid.TPC) continue;
      if (wids[w].Cryostat != tpcid.Cryostat)
        throw cet::exception("LineSegChanIntersect")
          << "Channel and line not in the same crostat.\n";

      // If the current wire id wire number is inbetween the start/end
      // point wires, the line segment intersects the wireID at some point.

      // TODO: for now, extend range, but that is application specific. fix asap
      // The longer we make the range, the more conservative it is, so it is safe
      // to extend the range a bit to get hits at the ends of the line
      unsigned int ext = 0;
      if (ExtendLine) ext = 10;

      if (lar::util::ValueInRange(wids[w].Wire * 1., (startW - ext) * 1., (endW + ext) * 1.))
        widsCrossed.push_back(wids[w]);
    }

    return widsCrossed.size() > 0;
  }

  //----------------------------------------------------------
  std::vector<double> APAGeometryAlg::ThreeChanPos(uint32_t u, uint32_t v, uint32_t z) const
  {

    // Say we've associated a U, V, and Z channel -- perhaps by associating hits
    // or cluster endpoints -- these don't necessarily all intersect, but they
    // are hopefully pretty close. Find the center of the 3 intersections.

    // get data needed along the way
    std::vector<geo::WireIDIntersection> UVIntersects;
    this->APAChannelsIntersect(u, v, UVIntersects);
    std::vector<double> UVzToZ(UVIntersects.size());
    geo::WireID Zwid = fGeom->ChannelToWire(z)[0];
    unsigned int cryo = Zwid.Cryostat;
    unsigned int tpc = Zwid.TPC;
    std::vector<geo::WireID> Uwids = fGeom->ChannelToWire(u);
    std::vector<geo::WireID> Vwids = fGeom->ChannelToWire(v);
    std::vector<geo::WireID> UwidsInTPC, VwidsInTPC;
    for (size_t i = 0; i < Uwids.size(); i++)
      if (Uwids[i].TPC == tpc) UwidsInTPC.push_back(Uwids[i]);
    for (size_t i = 0; i < Vwids.size(); i++)
      if (Vwids[i].TPC == tpc) VwidsInTPC.push_back(Vwids[i]);
    auto const Zcent = fGeom->WireIDToWireGeo(Zwid).GetCenter();

    std::cout << "Zcent = " << Zcent.Z() << ", UVintersects zpos = ";
    for (size_t uv = 0; uv < UVIntersects.size(); uv++) {
      std::cout << UVIntersects[uv].z << ", ";
    }
    std::cout << "\n";


    //std::cout << "U = " << u << " V = " << v << ", " << UVIntersects.size() << std::endl;

    if (UVIntersects.size() == 0) {
      if (UwidsInTPC.size() > 1 || VwidsInTPC.size() > 1)
        throw cet::exception("ThreeChanPos") << "U/V channels don't intersect, bad return.\n";

      // Now assume there are only one of each u and v wireIDs on in this TPC
      mf::LogWarning("ThreeChanPos")
        << "No U/V intersect, exceptional channels. See if U or V intersects Z\n";
      std::vector<double> yzCenter(2, 0.);
      geo::WireID Uwid = UwidsInTPC[0];
      geo::WireID Vwid = VwidsInTPC[0];
      geo::WireIDIntersection UZInt, VZInt;
      bool checkUZ = fGeom->WireIDsIntersect(Uwid, Zwid, UZInt);
      bool checkVZ = fGeom->WireIDsIntersect(Vwid, Zwid, VZInt);
      if (!checkUZ && !checkVZ)
        throw cet::exception("NoChanIntersect") << "No channels intersect, bad return.\n";
      if (checkUZ && !checkVZ) {
        yzCenter[0] = UZInt.y;
        yzCenter[1] = UZInt.z;
      }
      if (checkVZ && !checkUZ) {
        yzCenter[0] = VZInt.y;
        yzCenter[1] = VZInt.z;
      }
      if (checkUZ && checkVZ) {
        yzCenter[0] = (VZInt.y + UZInt.y) / 2;
        yzCenter[1] = (VZInt.z + UZInt.z) / 2;
      }
      return yzCenter;
    }


    // In case the uv channels intersect twice on the same side, choose the best case.
    // Note: this will not happen for APAs with UV angle at about 36, but will for 45
    std::cout << "UVzToZ = ";
    for (size_t widI = 0; widI < UVIntersects.size(); widI++) {
      UVzToZ[widI] = std::abs(UVIntersects[widI].z - Zcent.Z());
      std::cout << UVzToZ[widI] << ", ";
    }
    std::cout << "\n";

    unsigned int bestWidI = 0;
    double minZdiff = fGeom->TPC(Zwid).Length(); // start it out at maximum z
    for (unsigned int widI = 0; widI < UVIntersects.size(); widI++) {

      //std::cout << "widI = " << widI << std::endl;

      if (UVIntersects[widI].TPC == tpc && UVzToZ[widI] < minZdiff) {
        minZdiff = UVzToZ[widI];
        bestWidI = widI;
        //std::cout << "bestWidI = " << bestWidI << std::endl;
      }
    }
    geo::WireIDIntersection ChosenUVInt = UVIntersects[bestWidI];

    // Now having the UV intersection, get the UZ and VZ
    double UVInt[3] = {0.};
    UVInt[1] = ChosenUVInt.y;
    UVInt[2] = ChosenUVInt.z;
    geo::WireID Uwid = this->NearestWireIDOnChan(UVInt, u, 0, tpc, cryo);
    geo::WireID Vwid = this->NearestWireIDOnChan(UVInt, v, 1, tpc, cryo);
    geo::WireIDIntersection UZInt, VZInt;
    bool checkUZ = fGeom->WireIDsIntersect(Uwid, Zwid, UZInt);
    bool checkVZ = fGeom->WireIDsIntersect(Vwid, Zwid, VZInt);

    std::cout << "UZint.z = " << UZInt.z << " (" << checkUZ << "), VZint.z = " << VZInt.z << " ("
              << checkVZ << ")\n";

    // find the center
    std::vector<double> yzCenter(2, 0.);

    if (!checkUZ || !checkVZ) {
      //throw cet::exception("ThreeChanPos")  << "WireIDs were expected to intersect.\n";

      std::cout << "ChosenUVint.y = " << ChosenUVInt.y << "ChosenUVint.z = " << ChosenUVInt.z
                << std::endl;

      //temporary case
      yzCenter[0] = ChosenUVInt.y;
      yzCenter[1] = ChosenUVInt.z;
    }
    else {

      yzCenter[0] = (ChosenUVInt.y + UZInt.y + VZInt.y) / 3;
      yzCenter[1] = (ChosenUVInt.z + UZInt.z + VZInt.z) / 3;
    }

    return yzCenter;
  }

  //----------------------------------------------------------
  bool APAGeometryAlg::APAChannelsIntersect(
    uint32_t chan1,
    uint32_t chan2,
    std::vector<geo::WireIDIntersection>& IntersectVector) const
  {

    // Get the WireIDs and view for each channel, make sure views are different
    geo::WireIDIntersection widIntersect;
    std::vector<geo::WireID> wids1 = fGeom->ChannelToWire(chan1);
    std::vector<geo::WireID> wids2 = fGeom->ChannelToWire(chan2);
    geo::View_t view1 = fGeom->View(chan1);
    geo::View_t view2 = fGeom->View(chan2);
    if (view1 == view2) {
      mf::LogWarning("APAChannelsIntersect")
        << "Comparing two channels in the same view, return false";
      return false;
    }
    if (wids1[0].Cryostat != wids2[0].Cryostat ||
        this->ChannelToAPA(chan1) != this->ChannelToAPA(chan2)) {
      throw cet::exception("APAChannelsIntersect")
        << "Comparing two channels in in different APAs: "
        << "channel " << chan1 << " in Cryo " << wids1[0].Cryostat << ", APA "
        << this->ChannelToAPA(chan1) << ", and channel " << chan2 << " in Cryo "
        << wids2[0].Cryostat << ", APA " << this->ChannelToAPA(chan2) << "\n";
      return false;
    }

    // Loop through wids1 and see if wids2 has any intersecting wires,
    // given that the WireIDs are in the same TPC
    for (unsigned int i1 = 0; i1 < wids1.size(); i1++) {
      for (unsigned int i2 = 0; i2 < wids2.size(); i2++) {

        // make sure it is reasonable to intersect
        if (wids1[i1].Plane == wids2[i2].Plane ||
            wids1[i1].TPC != wids2[i2].TPC || // not reasonable for a *WireID*
            wids1[i1].Cryostat != wids2[i2].Cryostat)
          continue;

        //      std::cout << "Checking: \n WireID 1 = ("
        //                << wids1[i1].Cryostat << "," << wids1[i1].TPC << ","
        //                << wids1[i1].Plane    << "," << wids1[i1].Wire
        //                << ") \n WireID 2 = ("
        //                << wids2[i2].Cryostat << "," << wids2[i2].TPC << ","
        //                << wids2[i2].Plane    << "," << wids2[i2].Wire << ")" << std::endl;

        // Check if they even intersect; if they do, push back
        if (fGeom->WireIDsIntersect(wids1[i1], wids2[i2], widIntersect)) {

          //      std::cout << "we have an intersect" << std::endl;

          IntersectVector.push_back(widIntersect);
        }
      }
    }

    // Of all considered configurations, there are never more than
    // 4 intersections per channel pair
    if (IntersectVector.size() > 4) {
      mf::LogWarning("APAChannelsIntersect")
        << "Got " << IntersectVector.size() << " intersections for channels " << chan1 << " and "
        << chan2 << " - never expect more than 4, so far";
    }

    // With increasing IntersectVector index, the WireID
    // vector indices of the intersecting wireIDs increase.
    //   This matches the direction in which ChannelToWire
    //   builds its output WireID vector in the APA/35t Alg
    std::sort(IntersectVector.begin(), IntersectVector.end());

    // return true if any intersection points were found
    if (IntersectVector.size() == 0)
      return false;
    else
      return true;
  }

} //end namespace apa