View.cxx

Go to the documentation of this file.

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

#include "DiscreteTracker/View.h"

#include "ChannelInfo/BadChanList.h"
#include "Geometry/Geometry.h"
#include "RecoBase/CellHit.h"

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

namespace dt
{

  //......................................................................
  View::View(art::PtrVector<rb::CellHit> chits, bool considerBadChans)
  {
    assert(!chits.empty());

    fGeoView = chits[0]->View();

    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<chaninfo::BadChanList> bcl;

    fChans.resize(geom->NPlanes());

    fPlanes = geom->GetPlanesByView(fGeoView);

    rb::Cluster clust(chits); // Just to get the extents
    // Clumsy technique to remove all the planes outside the hit range
    const unsigned int minPlane = clust.MinPlane();
    for(unsigned int p = 0; p < minPlane; ++p) fPlanes.erase(p);
    const unsigned int maxPlane = clust.MaxPlane();
    for(unsigned int p = maxPlane+1; p < geom->NPlanes(); ++p)
      fPlanes.erase(p);

    // Have the box we save be one bigger than the extent of the hit cells (but
    // not larger than the detector).
    const unsigned int minCell = std::max(clust.MinCell(fGeoView), 1u)-1;
    fMinCell = minCell;
    const unsigned int maxCell = std::min(clust.MaxCell(fGeoView)+1, geom->Plane(minPlane)->Ncells()-1);

    for(std::set<unsigned int>::iterator it = fPlanes.begin();
        it != fPlanes.end(); ++it){
      const unsigned int plane = *it;

      // Are you running with a geometry inconsistent with what you generated
      // with?
      assert(plane < geom->NPlanes());

      fChans[plane].reserve(maxCell-minCell+1);

      for(unsigned int cell = minCell; cell <= maxCell; ++cell){
        dt::Channel chan(plane, cell);
        if(considerBadChans &&
           (bcl->IsBad(plane, cell) || bcl->IsLowEfficiency(plane, cell))){
          chan.SetDead();
        }

        fChans[plane].push_back(chan);
      }
    } // end for it (plane)

    for(unsigned int n = 0; n < chits.size(); ++n){
      const art::Ptr<rb::CellHit> chit = chits[n];
      const int plane = chit->Plane();
      const int cell = chit->Cell();
      dt::Channel& chan = fChans[plane][cell-minCell];
      // Behave correctly in the face of eg weird CalHit configuration, and of
      // channels that have bad efficiency but somehow worked this one time.
      if(chan.Type() != Channel::kDead ||
         (considerBadChans && bcl->IsLowEfficiency(plane, cell))){
        chan.SetHit(chit);
      }
    }

    // Don't consider any planes that are completely dead
    std::vector<unsigned int> deadPlanes;
    for(std::set<unsigned int>::iterator it = fPlanes.begin();
        it != fPlanes.end(); ++it){
      const unsigned int plane = *it;

      bool alldead = true;
      for(unsigned int cell = 0; cell < fChans[plane].size(); ++cell){
        if(fChans[plane][cell].Type() != Channel::kDead){
          alldead = false;
          break;
        }
      } // end for cell
      if(alldead) deadPlanes.push_back(plane);
    } // end for it (plane)

    for(unsigned int n = 0; n < deadPlanes.size(); ++n)
      fPlanes.erase(deadPlanes[n]);
  }

  //......................................................................
  int View::PrevPlane(unsigned int plane) const
  {
    std::set<unsigned int>::iterator it = fPlanes.find(plane);
    if(it == fPlanes.end()) return -1;
    if(it == fPlanes.begin()) return -1;
    --it;
    return *it;
  }

  //......................................................................
  int View::NextPlane(unsigned int plane) const
  {
    std::set<unsigned int>::iterator it = fPlanes.find(plane);
    if(it == fPlanes.end()) return -1;
    ++it;
    if(it == fPlanes.end()) return -1;
    return *it;
  }

  //......................................................................
  int View::AdjacentPlane(unsigned int plane, Direction dir) const
  {
    if(dir == kUpstream) return PrevPlane(plane);
    return NextPlane(plane);
  }

  //......................................................................
  dt::Channel View::Channel(int plane, int cell) const
  {
    if(Outside(plane, cell)){
      dt::Channel outside(plane, cell);
      outside.SetDead();
      return outside;
    }

    return fChans[plane][cell-fMinCell];
  }

  //......................................................................
  dt::Channel& View::Channel(int plane, int cell)
  {
    assert(!Outside(plane, cell));

    return fChans[plane][cell-fMinCell];
  }

  //......................................................................
  View::ChunkMap View::MakeChunks() const
  {
    ChunkMap chunks_map;

    for(std::set<unsigned int>::iterator it = fPlanes.begin();
        it != fPlanes.end(); ++it){
      const int plane = *it;

      bool inside = false;
      int begin = -1;

      const unsigned int cellMax = fMinCell+fChans[plane].size();
      for(unsigned int cell = fMinCell; cell < cellMax; ++cell){
        const Channel::Type_t type = Channel(plane, cell).Type();
        if(!inside && (type == Channel::kHit || type == Channel::kDead)){
          inside = true;
          begin = cell;
        }
        if(inside && (type != Channel::kHit && type != Channel::kDead)){
          inside = false;
          std::vector<Chunk> chunks = MakeChunkCombos(plane, begin, cell);
          chunks_map[plane].insert(chunks_map[plane].begin(),
                                   chunks.begin(), chunks.end());
        }
      } // end for cell
      if(inside){
        std::vector<Chunk> chunks = MakeChunkCombos(plane, begin, cellMax);
        chunks_map[plane].insert(chunks_map[plane].begin(),
                                 chunks.begin(), chunks.end());
      }
    } // end for it (plane)

    return chunks_map;
  }

  //......................................................................
  std::vector<Chunk> View::MakeChunkCombos(int plane,
                                           int begin, int end) const
  {
    std::vector<Chunk> ret;

    assert(begin < end);

    // Pairs of hit cell and extremal neighbouring dead cell, ordered low to high
    std::vector<std::pair<int, int> > lows, highs;

    // Every hit cell that has a boundary below, and every one with a boundary above
    for(int k = begin; k < end; ++k){
      const Channel::Type_t here = Channel(plane, k).Type();
      const Channel::Type_t below = Channel(plane, k-1).Type();
      const Channel::Type_t above = Channel(plane, k+1).Type();

      if(here == Channel::kHit && (k == begin || below == Channel::kDead)){
        int lower = k;
        // Run the dead cell along until we find the end of the dead region
        while(Channel(plane, lower-1).Type() == Channel::kDead && lower-1 >= begin) --lower;
        lows.push_back(std::make_pair(lower, k));
      }

      if(here == Channel::kHit && (k == end-1 || above == Channel::kDead)){
        int higher = k;
        while(Channel(plane, higher+1).Type() == Channel::kDead && higher+1 < end) ++higher;
        highs.push_back(std::make_pair(k, higher));
      }
    } // end for k


    // Every combination of those starting and stopping points
    for(unsigned int i = 0; i < lows.size(); ++i){
      for(unsigned int j = 0; j < highs.size(); ++j){
        if(lows[i].second > highs[j].first) continue;

        art::PtrVector<rb::CellHit> hits;
        for(int k = lows[i].second; k <= highs[j].first; ++k){
          if(Channel(plane, k).Type() == Channel::kHit)
            hits.push_back(Channel(plane, k).GetHit());
        }

        for(bool up: {false, true}){
          Chunk chunk;
          if(up) chunk = Chunk(plane, lows[i].first, lows[i].second, highs[j].first, highs[j].second, up, hits);
          else   chunk = Chunk(plane, highs[j].first, highs[j].second, lows[i].first, lows[i].second, up, hits);

          ret.push_back(chunk);
        } // end for up
      } // end for i
    } // end for j


    // Also, make a chunk for every contiguous block of dead channels. A track
    // can go through these without needing to hit the live cells
    // adjacent. Hopefully the scoring system weighting real hits higher means
    // this won't cause problems.
    for(int i = begin; i < end; ++i){
      if(Channel(plane, i).Type() != Channel::kDead) continue;
      const int from = i;
      for(; i < end; ++i){
        if(Channel(plane, i).Type() != Channel::kDead) break;
      }
      const int to = i-1;

      for(bool up: {false, true}){
        const Chunk chunk(plane, from, to, from, to, up,
                          art::PtrVector<rb::CellHit>());
        ret.push_back(chunk);
      }
    }

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

  //......................................................................
  bool View::Outside(int plane, int cell) const
  {
    return (plane < 0 || cell < int(fMinCell) ||
            plane >= int(fChans.size()) ||
            cell-int(fMinCell) >= int(fChans[plane].size()));
  }

} // namespace