Cand.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    Cand.cxx
// \version $Id: Cand.cxx,v 1.28 2012-10-31 00:05:29 bckhouse Exp $
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "DiscreteTracker/Cand.h"

#include "DiscreteTracker/Chain.h"

#include <algorithm>
#include <cassert>

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"

#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"

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

namespace dt
{
  //......................................................................
  Cand::Cand(geo::View_t view, bool up)
    : fView(view), fUp(up),
      fAllMMax(+1e10), fAllMMin(-1e10),
      fAllCMax(+1e10), fAllCMin(-1e10),
      fNChunks(0), fNDeadChunks(0),
      fAllHitsDirty(true)
  {
  }

  //......................................................................
  void Cand::AddChunk(const Chunk& chunk, Direction dir)
  {
    fAllHitsDirty = true;

    if(chunk.HasHits())
      ++fNChunks;
    else
      ++fNDeadChunks;
    assert(chunk.Up() == fUp);

    if(dir == kDownstream)
      fChunks.push_back(chunk);
    else
      fChunks.push_front(chunk);
  }

  //......................................................................
  void Cand::PopFirstChunk()
  {

    if(fChunks.front().HasHits())
      --fNChunks;
    else
      --fNDeadChunks;

    fChunks.pop_front();

    fAllHitsDirty = true;
  }

  //......................................................................
  void Cand::PopLastChunk()
  {
    if(fChunks.back().HasHits())
      --fNChunks;
    else
      --fNDeadChunks;

    fChunks.pop_back();

    fAllHitsDirty = true;
  }

  //......................................................................
  bool Cand::operator<(const Cand& rhs) const
  {
    // This is sufficient to sort them, but can't be used for a map or set
    // Possibly shouldn't include it as an operator then...

    if(NHitPlanes() < rhs.NHitPlanes()) return true;
    if(NHitPlanes() > rhs.NHitPlanes()) return false;

    if(fSegs.size() < rhs.fSegs.size()) return true;
    if(fSegs.size() > rhs.fSegs.size()) return false;

    // Having more dead is a _bad_ thing. Select the ones without
    // preferentially
    if(fNDeadChunks > rhs.fNDeadChunks) return true;
    if(fNDeadChunks < rhs.fNDeadChunks) return false;

    //    if(fChunks.size() < rhs.fChunks.size()) return true;
    //    if(fChunks.size() > rhs.fChunks.size()) return false;

    // Even if these are only destined to be singleton chunks, need to give the
    // biggest ones back first because we'll stop looking for new chains when
    // we see small ones.
    if(AllHits().size() < rhs.AllHits().size()) return true;
    if(AllHits().size() > rhs.AllHits().size()) return false;

    // Arbitrary, but removes dependence on order FindCands emits things
    return FirstChunk().Plane() < rhs.FirstChunk().Plane();
  }

  //......................................................................
  void Cand::MinMaxLines(const SubSeg& segi, const SubSeg& segj,
                         double& mmax, double& mmin,
                         double& cmax, double& cmin) const
  {
    const double dz = segj.z-segi.z;

    // Got sorted at some earlier stage
    assert(dz > 0);

    mmax = std::min(mmax, (segj.y1-segi.y0)/dz);
    mmin = std::max(mmin, (segj.y0-segi.y1)/dz);

    cmax = std::min(cmax, segi.y1-segi.z*(segj.y0-segi.y1)/dz);
    cmin = std::max(cmin, segi.y0-segi.z*(segj.y1-segi.y0)/dz);
  }

  //......................................................................
  void Cand::MinMaxLines(const Segment& segi, const Segment& segj,
                         double& mmax, double& mmin,
                         double& cmax, double& cmin) const
  {
    const double dz = segj.zL-segi.zL;

    if(dz < 0){
      MinMaxLines(segj, segi, mmax, mmin, cmax, cmin);
      return;
    }

    mmax = +1e10;
    mmin = -1e10;
    cmax = +1e10;
    cmin = -1e10;

    MinMaxLines(segi.GetSubSeg(kUpstream),   segj.GetSubSeg(kUpstream),   mmax, mmin, cmax, cmin);
    MinMaxLines(segi.GetSubSeg(kUpstream),   segj.GetSubSeg(kDownstream), mmax, mmin, cmax, cmin);
    MinMaxLines(segi.GetSubSeg(kDownstream), segj.GetSubSeg(kUpstream),   mmax, mmin, cmax, cmin);
    MinMaxLines(segi.GetSubSeg(kDownstream), segj.GetSubSeg(kDownstream), mmax, mmin, cmax, cmin);
  }

  //......................................................................
  bool Cand::TryAddSeg(const Segment& seg, Direction dir, bool check)
  {
    if(check){
      // Must alternate sides
      assert(fSegs.empty() || seg.left != ExtremalSeg(dir).left);
    }

    double newAllMMax = fAllMMax;
    double newAllMMin = fAllMMin;
    double newAllCMax = fAllCMax;
    double newAllCMin = fAllCMin;

    for(std::list<Segment>::iterator it = fSegs.begin(); it != fSegs.end(); ++it){
      const Segment& segi = *it;
      // In principle would have to loop over all j > i. But this happens every
      // time we add a new seg, so only the last one is novel.
      const Segment& segj = seg;

      double mmax, mmin, cmax, cmin;
      MinMaxLines(segi, segj, mmax, mmin, cmax, cmin);
      newAllMMax = std::min(mmax, newAllMMax);
      newAllMMin = std::max(mmin, newAllMMin);
      newAllCMax = std::min(cmax, newAllCMax);
      newAllCMin = std::max(cmin, newAllCMin);

      if(newAllMMax <= newAllMMin ||
         newAllCMax <= newAllCMin) return false;
    } // end for it

    // TODO: I'm worried about the interaction of this with imperfect
    // geometry/staggering, but it, or something like it, is necessary to
    // prevent cands with the obviously "wrong" sense through one-cell chunks
    if((fUp && newAllMMax < 0) || (!fUp && newAllMMin > 0)) return false;

    // Everything checks out, finalize the addition

    if(dir == kDownstream)
      fSegs.push_back(seg);
    else
      fSegs.push_front(seg);

    fAllHitsDirty = true;

    fAllMMax = newAllMMax;
    fAllMMin = newAllMMin;
    fAllCMax = newAllCMax;
    fAllCMin = newAllCMin;

    return true;
  }

  //......................................................................
  bool Cand::TryAddSegOutOfOrder(const Segment& seg)
  {
    if(TryAddSeg(seg, kDownstream, false)){
      fSegs.sort();
      return true;
    }
    return false;
  }

  //......................................................................
  void Cand::PopFirstSeg()
  {
    assert(fSegs.size() > 1);

    const Segment segi = fSegs.front();
    fSegs.pop_front();

    PopSegHelper(segi, kUpstream);
  }

  //......................................................................
  void Cand::PopLastSeg()
  {
    assert(fSegs.size() > 1);

    const Segment segi = fSegs.back();
    fSegs.pop_back();

    PopSegHelper(segi, kDownstream);
  }

  //......................................................................
  void Cand::PopSegHelper(const Segment& segi, Direction end)
  {
    if(fSegs.size() == 1){
      fAllMMax = +1e10;
      fAllMMin = -1e10;
      fAllCMax = +1e10;
      fAllCMin = -1e10;
      return;
    }

    // Work out if this seg was responsible for any of the strongest
    // constraints.
    bool needRecalc = false;
    for(std::list<Segment>::iterator it = fSegs.begin(); it != fSegs.end(); ++it){
      const Segment& segj = *it;

      double mmax, mmin, cmax, cmin;
      MinMaxLines((end == kDownstream) ? segj : segi,
                  (end == kDownstream) ? segi : segj,
                  mmax, mmin, cmax, cmin);

      if(mmax-1e-6 < fAllMMax ||
         mmin+1e-6 > fAllMMin ||
         cmax-1e-6 < fAllCMax ||
         cmin+1e-6 > fAllCMin){
        needRecalc = true;
        break;
      }
    }

    if(!needRecalc) return;

    // If this seg was important, need to recalculate the whole lot to find the
    // 2nd most stringent. To avoid this we'd have to store more information in
    // the object.
    std::list<Segment> segs = fSegs;
    fSegs.clear();

    fAllMMax = +1e10;
    fAllMMin = -1e10;
    fAllCMax = +1e10;
    fAllCMin = -1e10;

    for(std::list<Segment>::iterator it = segs.begin(); it != segs.end(); ++it){
      const bool ok = TryAddSeg(*it, kDownstream);
      assert(ok);
    }
  }

  //......................................................................
  void Cand::EstimateStraightLine(double& z1, double& y1,
                                  double& z2, double& y2) const
  {
    assert(fSegs.size() > 1);

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

    // Best guess is starting halfway to the previous plane and ending halfway
    // to the next.
    const int firstPlane = FirstChunk().Plane();
    const int prevPlane = fSegs.front().left ? geom->NextPlaneInView(firstPlane, -1) : firstPlane;

    const int lastPlane = LastChunk().Plane();
    const int nextPlane = fSegs.back().left ? lastPlane : geom->NextPlaneInView(lastPlane, +1);

    if(prevPlane != geo::kPLANE_NOT_FOUND){
      double a[3], b[3];
      geom->Plane(prevPlane)->Cell(0)->GetCenter(a);
      geom->Plane(firstPlane)->Cell(0)->GetCenter(b);
      z1 = (a[2]+b[2])/2;
    }
    else{
      z1 = fSegs.front().zL;
    }

    if(nextPlane != geo::kPLANE_NOT_FOUND){
      double a[3], b[3];
      geom->Plane(lastPlane)->Cell(0)->GetCenter(a);
      geom->Plane(nextPlane)->Cell(0)->GetCenter(b);
      z2 = (a[2]+b[2])/2;
    }
    else{
      z2 = fSegs.back().zR;
    }

    // Except if the first or last chunk is partial, in which case our best
    // guess is nearest the centre of the first or last cell.
    if(!fSegs.front().left){
      const double newz1 = ClosestToEndCell(true);
      if(newz1 < fSegs.front().zL) z1 = newz1;
    } // end if


    if(fSegs.back().left){
      const double newz2 = ClosestToEndCell(false);
      if(newz2 > fSegs.back().zR) z2 = newz2;
    }

    const double cavg = CAvg();
    const double mavg = MAvg();

    y1 = cavg+mavg*z1;
    y2 = cavg+mavg*z2;

    assert(z2 >= z1);
  }

  //......................................................................
  int Cand::NHitPlanes() const
  {
    int ret = 0;
    for(std::list<Chunk>::const_iterator it = fChunks.begin(); it != fChunks.end(); ++it){
      if(!it->AllDead()) ++ret;
    }

    return ret;
  }

  //......................................................................
  const Chunk& Cand::FirstChunk() const
  {
    return fChunks.front();
  }

  //......................................................................
  const Chunk& Cand::LastChunk() const
  {
    return fChunks.back();
  }

  //......................................................................
  const Segment& Cand::FirstSeg() const
  {
    return fSegs.front();
  }

  //......................................................................
  const Segment& Cand::LastSeg() const
  {
    return fSegs.back();
  }

  //......................................................................
  bool Cand::IsExtremalChunkComplete(Direction dir) const
  {
    if(dir == kUpstream) return FirstSeg().left;
    return !LastSeg().left;
  }

  //......................................................................
  void Cand::MarkExtremalChunkShower(Direction dir)
  {
    const bool complete = IsExtremalChunkComplete(dir);

    if(dir == kUpstream){
      fChunks.begin()->SetShowerChunk();
      if(complete) PopFirstSeg();
      PopFirstSeg();
      const std::pair<Segment, Segment> segs = fChunks.begin()->GetSegs();
      bool ok = TryAddSeg(segs.second, kUpstream);
      assert(ok);
      if(complete){
        ok = TryAddSeg(segs.first, kUpstream);
        assert(ok);
      }
    }
    else{
      fChunks.rbegin()->SetShowerChunk();
      if(complete) PopLastSeg();
      PopLastSeg();
      const std::pair<Segment, Segment> segs = fChunks.rbegin()->GetSegs();
      bool ok = TryAddSeg(segs.first, kDownstream);
      assert(ok);
      if(complete){
        ok = TryAddSeg(segs.second, kDownstream);
        assert(ok);
      }
    }
  }

  //......................................................................
  bool Cand::IsAllShower() const
  {
    for(std::list<Chunk>::const_iterator it = fChunks.begin(); it != fChunks.end(); ++it){
      if(!it->IsShowerChunk()) return false;
    }
    return true;
  }

  //......................................................................
  art::PtrVector<rb::CellHit> Cand::AllHits() const
  {
    if(!fAllHitsDirty) return fAllHits;
    fAllHitsDirty = false;

    fAllHits.clear();
    if(fSegs.size() <= 1) return fAllHits;
    if(fNChunks == 0) return fAllHits;

    for(std::list<Chunk>::const_iterator it = fChunks.begin(); it != fChunks.end(); ++it){
      art::PtrVector<rb::CellHit> chunkHits = it->AllHits();

      // Unless it's the type of chunk you have to do by dead-reckoning the
      // track.
      if(it->IsShowerChunk() ||
         (it == fChunks.begin() && !FirstSeg().left) ||
         (it->Plane() == LastChunk().Plane() && LastSeg().left)){

        chunkHits = it->HitsOnLine(MAvg(), CAvg());
      }

      for(unsigned int m = 0; m < chunkHits.size(); ++m)
        fAllHits.push_back(chunkHits[m]);
    } // end for it

    return fAllHits;
  }

  //......................................................................
  rb::Track Cand::ToTrack() const
  {
    double z1, y1, z2, y2;
    EstimateStraightLine(z1, y1, z2, y2);
    rb::Track ret(art::PtrVector<rb::CellHit>(),
                  fView, y1, z1, y2-y1, z2-z1);
    ret.AppendTrajectoryPoint(y2, z2);
    for(std::list<Chunk>::const_iterator it = fChunks.begin(); it != fChunks.end(); ++it){
      ret.Add(it->AllHits());
    }
    // TODO: this can happen sometimes, figure out when and why
    //    assert(ret.NCell() > 0);
    return ret;
  }

  //......................................................................
  double Cand::ClosestToEndCell(bool left) const
  {
    double bestL = left ? +1e10 : -1e10;
    double bestZ = -1e10;

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

    const art::PtrVector<rb::CellHit> allhits = AllHits();
    for(unsigned int n = 0; n < allhits.size(); ++n){
      const geo::CellGeo* geocell = geom->Plane(allhits[n]->Plane())->Cell(allhits[n]->Cell());
      double xyz[3];
      geocell->GetCenter(xyz);
      // y-view by convention
      xyz[1] = xyz[fView];
      xyz[0] = 0;

      const TVector3 orig(0, CAvg(), 0);
      const TVector3 dir(0, MAvg(), 1);

      TVector3 start;
      geo::ClosestApproach(TVector3(xyz), orig, dir, start);

      // Distance along the line
      const double L = (start-orig).Dot(dir);

      if((left && L < bestL) ||
         (!left && L > bestL)){
        bestL = L;
        bestZ = start.Z();
      }
    } // end for n

    // TODO
    //    assert(bestZ > -1e10);
    return bestZ;
  }

  //......................................................................
  void Cand::TrimFront()
  {
    fAllHitsDirty = true;

    while(!fChunks.empty()){
      if(!fChunks.front().AllDead()) return;
      const int plane = fChunks.front().Plane();
      fChunks.pop_front();
      while(fSegs.front().plane == plane) fSegs.pop_front();
    }
  }

  //......................................................................
  void Cand::TrimBack()
  {
    fAllHitsDirty = true;

    while(!fChunks.empty()){
      if(!fChunks.back().AllDead()) return;
      const int plane = fChunks.back().Plane();
      fChunks.pop_back();
      while(fSegs.back().plane == plane) fSegs.pop_back();
    }
  }

  //......................................................................
  void Cand::ShortenOne(Direction dir)
  {
    if(dir == kDownstream){
      if(!FirstSeg().left) PopFirstChunk();
      PopFirstSeg();
    }
    else{
      if(LastSeg().left) PopLastChunk();
      PopLastSeg();
    }
  }


  //......................................................................
  Cand Cand::MaybeFlip(bool& ok) const
  {
    Cand ret(fView, !fUp);

    for(std::list<Chunk>::const_iterator it = fChunks.begin(); it != fChunks.end(); ++it){
      Chunk chunk = *it;
      chunk.SetUp(!fUp);
      ret.AddChunk(chunk, kDownstream);

      const std::pair<Segment, Segment> segs = chunk.GetSegs();

      // Add first seg unless it's the very first and we don't have it
      if(it != fChunks.begin() || fSegs.front().left){
        ok = ret.TryAddSeg(segs.first, kDownstream);
        if(!ok) return ret;
      }

      // Add second seg unless it's the last and we don't have it
      if(ret.fSegs.size() < fSegs.size()){
        ok = ret.TryAddSeg(segs.second, kDownstream);
        if(!ok) return ret;
      }
    }

    return ret;
  }

  //......................................................................
  bool Cand::FindChunk(int plane, Chunk& chunk) const
  {
    for(std::list<Chunk>::const_iterator it = fChunks.begin(); it != fChunks.end(); ++it){
      if(it->Plane() == plane){
        chunk = *it;
        return true;
      }
    }
    return false;
  }

  //......................................................................
  bool Cand::operator==(const Cand& rhs) const
  {
    if(fView != rhs.fView) return false;
    if(fUp != rhs.fUp) return false;

    return fSegs == rhs.fSegs;
  }

} // namespace