Chain.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    Chain.cxx
// \version $Id: Chain.cxx,v 1.27 2012-11-01 06:05:03 bckhouse Exp $
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "DiscreteTracker/Chain.h"

#include "DiscreteTracker/Cand.h"

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

#include "RecoBase/CellHit.h"
#include "RecoBase/HitMap.h"
#include "RecoBase/Track.h"

#include "Utilities/func/MathUtil.h"

namespace dt
{
  //......................................................................
  Chain::Chain(const Cand& c)
  {
    Add(c);
  }

  //......................................................................
  rb::Track Chain::ToTrack() const
  {
    art::ServiceHandle<geo::Geometry> geom;

    // Special case for vertical cands (who'll otherwise estimate themselves
    // arbitrarily at a slope depending if they're marked "up" or "down").
    if(ExtentPlane() == 0 && AllHits().size() > 1){
      const Chunk& chunk = FirstCand().FirstChunk();
      double xyz0[3];
      geom->Plane(chunk.Plane())->Cell(chunk.TopHitCell())->GetCenter(xyz0);
      double xyz1[3];
      geom->Plane(chunk.Plane())->Cell(chunk.BottomHitCell())->GetCenter(xyz1);

      geo::View_t view = FirstCand().View();
      rb::Track tr(AllHits(), FirstCand().View(),
                   xyz0[view], xyz0[2],
                   xyz1[view]-xyz0[view], xyz1[2]-xyz0[2]);

      tr.AppendTrajectoryPoint(xyz1[view], xyz1[2]);
      return tr;
    }

    // The algorithm is to take the weighted average of the y-positions of all
    // relevant candidates at each z-position. (Convention is that the non-z
    // coordinate is called "y").

    // Estimate the start and end points of all the component cands upfront
    std::vector<TVector3> p0s, p1s;
    for(std::list<Cand>::const_iterator candIt = fCands.begin();
        candIt != fCands.end(); ++candIt){
      double z1, y1, z2, y2;
      candIt->EstimateStraightLine(z1, y1, z2, y2);
      p0s.push_back(TVector3(0, y1, z1));
      p1s.push_back(TVector3(0, y2, z2));
    }

    std::vector<TVector3> pts;

    // TODO, figure out where -1e10 Z's come from
    const double minZ = std::max(0., p0s[0].Z());
    const double maxZ = std::min(p1s.back().Z(), geom->DetLength());

    const double kZStep = 5; // cm
    const int stepMax = int((maxZ-minZ)/kZStep);

    for(int step = 0; step <= stepMax; ++step){
      const double z = minZ+(double(step)/stepMax)*(maxZ-minZ);

      TVector3 meanp;
      double totweight = 0;

      for(unsigned int pIdx = 0; pIdx < p0s.size(); ++pIdx){
        const TVector3& p0 = p0s[pIdx];
        const TVector3& p1 = p1s[pIdx];
        const double dz = p1.Z()-p0.Z();

        const TVector3 p = p0*((p1.Z()-z)/dz) + p1*((z-p0.Z())/dz);

        // Give more weight to cands near their centres
        double weight = (z-p0.Z())/dz;
        if(weight > 0.5) weight = 1-weight;
        weight += 1e-3; // Avoid zeros at very ends
        if(weight < 0) continue;

        meanp += p*weight;
        totweight += weight;
      }

      // This can happen when two cands are connected because they intersect at
      // a shallow angle inside a large dead region. There's probably a neater
      // way to illustrate this in the trajectory, but it's rare, so this will
      // do.
      if(totweight == 0) continue; // No cands cover this range

      meanp *= 1/totweight;

      pts.push_back(meanp);
    } // end for step

    // TODO
    //    assert(pts.size() >= 2);
    assert(!pts.empty());
    if(pts.size() == 1){
      mf::LogWarning("DiscreteTracker") << "Only one trajectory point. Fabricating a second one." << std::endl;
      pts.push_back(pts[0]+TVector3(0, 0, 1));
    }

    const TVector3 dir = pts[1]-pts[0];

    rb::Track tr(AllHits(),
                 fCands.begin()->View(),
                 pts[0].Y(), pts[0].Z(), dir.Y(), dir.Z());

    for(unsigned int n = 1; n < pts.size(); ++n)
      tr.AppendTrajectoryPoint(pts[n].Y(), pts[n].Z());

    return tr;
  }

  //......................................................................
  std::vector<rb::Track> Chain::ToDebugTracks() const
  {
    std::vector<rb::Track> ret;

    for(std::list<Cand>::const_iterator candIt = fCands.begin();
        candIt != fCands.end(); ++candIt){
      ret.push_back(candIt->ToTrack());
    }

    return ret;
  }

  //......................................................................
  void Chain::Add(const Cand& c)
  {
    assert(c.NChunks() > 0);
    assert(c.NSegs() > 1);
    fCands.push_back(c);
  }

  //......................................................................
  void Chain::Add(const Cand& c, Direction dir)
  {
    if(dir == kDownstream)
      Add(c);
    else
      AddFront(c);
  }

  //......................................................................
  void Chain::AddFront(const Cand& c)
  {
    assert(c.NChunks() > 0);
    assert(c.NSegs() > 1);
    fCands.push_front(c);
  }

  //......................................................................
  unsigned int Chain::FirstPlane() const
  {
    return FirstCand().FirstChunk().Plane();
  }

  //......................................................................
  unsigned int Chain::LastPlane() const
  {
    return LastCand().LastChunk().Plane();
  }

  //......................................................................
  unsigned int Chain::ExtremalPlane(Direction dir) const
  {
    return ExtremalCand(dir).ExtremalPlane(dir);
  }

  //......................................................................
  unsigned int Chain::ExtremalHitPlane(Direction dir) const
  {
    std::vector<Chunk> chunks = AllChunks();
    if(dir == kUpstream) std::reverse(chunks.begin(), chunks.end());

    while(!chunks.empty() && chunks.back().AllDead()) chunks.pop_back();

    return chunks.back().Plane();
  }

  //......................................................................
  const Chunk& Chain::FirstChunk() const
  {
    return FirstCand().FirstChunk();
  }

  //......................................................................
  const Chunk& Chain::LastChunk() const
  {
    return LastCand().LastChunk();
  }

  //......................................................................
  const Segment& Chain::FirstSeg() const
  {
    return FirstCand().FirstSeg();
  }

  //......................................................................
  const Segment& Chain::LastSeg() const
  {
    return LastCand().LastSeg();
  }

  //......................................................................
  art::PtrVector<rb::CellHit> Chain::AllHits() const
  {
    art::PtrVector<rb::CellHit> ret;

    // To eliminate duplicates
    rb::HitMap hmap;

    for(std::list<Cand>::const_iterator it = fCands.begin();
        it != fCands.end(); ++it){
      const art::PtrVector<rb::CellHit>& hits = it->AllHits();
      for(unsigned int hitIdx = 0; hitIdx < hits.size(); ++hitIdx){
        const art::Ptr<rb::CellHit> hit = hits[hitIdx];
        if(!hmap.CellExists(hit->Plane(), hit->Cell())){
          ret.push_back(hit);
        }
        hmap.Add(hit);
      }
    } // end for it

    // TODO
    //    assert(!ret.empty());
    return ret;
  }

  //......................................................................
  std::vector<Chunk> Chain::AllChunks() const
  {
    // TODO: this doesn't take into account the possibility that chunks on the
    // same plane don't entirely correspond (due to dead channels?)
    std::map<unsigned int, Chunk> retMap;

    for(std::list<Cand>::const_iterator it = fCands.begin(); it != fCands.end(); ++it){
      const std::list<Chunk> chunks = it->AllChunks();
      for(std::list<Chunk>::const_iterator it = chunks.begin(); it != chunks.end(); ++it){
        retMap[it->Plane()] = *it;
      }
    }

    // Come out sorted in plane order
    std::vector<Chunk> ret;
    for(std::map<unsigned int, Chunk>::iterator it = retMap.begin(); it != retMap.end(); ++it)
      ret.push_back(it->second);
    return ret;
  }

  //......................................................................
  int Chain::NHitPlanes() const
  {
    std::vector<Chunk> chunks = AllChunks();
    int ret = 0;
    for(unsigned int n = 0; n < chunks.size(); ++n)
      if(!chunks[n].AllDead()) ++ret;

    return ret;
  }

  //......................................................................
  void Chain::TrimEnds()
  {
    const unsigned int origSize = AllHits().size();

    assert(origSize > 0);

    while(true){
      const Cand cand = fCands.front();
      fCands.pop_front();
      if(fCands.empty() || AllHits().size() < origSize){
        fCands.push_front(cand);
        break;
      }
    }

    while(true){
      const Cand cand = fCands.back();
      fCands.pop_back();
      if(fCands.empty() || AllHits().size() < origSize){
        fCands.push_back(cand);
        break;
      }
    }

    fCands.front().TrimFront();
    fCands.back().TrimBack();
  }

  //......................................................................
  void Chain::Truncate(unsigned int plane, Direction end)
  {
    if(end == kUpstream){
      while(FirstPlane() < plane){
        fCands.begin()->PopFirstChunk();
        if(FirstCand().AllChunks().empty()) fCands.pop_front();
      }
      while(FirstSeg().plane < int(plane)){
        if(fCands.begin()->NSegs() > 1)
          fCands.begin()->PopFirstSeg();
        else
          fCands.pop_front();
      }
    }
    else{
      while(LastPlane() > plane){
        fCands.rbegin()->PopLastChunk();
        if(LastCand().AllChunks().empty()) fCands.pop_back();
      }
      while(LastSeg().plane > int(plane)){
        if(fCands.rbegin()->NSegs() > 1)
          fCands.rbegin()->PopLastSeg();
        else
          fCands.pop_back();
      }
    }

    if(!AllHits().empty()) TrimEnds();
  }

  //......................................................................
  dt::View::ChunkMap Chain::AsChunkMap() const
  {
    dt::View::ChunkMap ret;

    const std::vector<Chunk> chunks = AllChunks();
    for(unsigned int n = 0; n < chunks.size(); ++n)
      ret[chunks[n].Plane()].push_back(chunks[n]);

    return ret;
  }

  //......................................................................
  bool compareByLength(const Chain& a, const Chain& b)
  {
    return a.ExtentPlane() < b.ExtentPlane();
  }

  //......................................................................
  bool compareByStart(const Chain& a, const Chain& b)
  {
    return a.FirstPlane() < b.FirstPlane();
  }

  //......................................................................
  bool compareByEnd(const Chain& a, const Chain& b)
  {
    return a.LastPlane() < b.LastPlane();
  }

} // namespace