ElasticArmsHS_module.cc

Go to the documentation of this file.

///
/// \file    ElasticArmsHS.cxx
/// \brief   Application of the elastic arms fitter to hough results
/// \author  messier@indiana.edu
/// \version $Id: ElasticArmsHS.cxx,v 1.21 2012-10-09 01:15:52 edniner Exp $
///
#include "ElasticArms/func/ElasticArms.h"
#include "ElasticArms/func/Minimizer.h"
#include "ElasticArms/art/GridSearch.h"

#include <algorithm>

#include "TNtuple.h"
#include "TH1F.h"

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "art/Framework/Services/Optional/TFileService.h"

#include "ChannelInfo/BadChanList.h"
#include "GeometryObjects/CellGeo.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/PlaneGeo.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Prong.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/HoughResult.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"

/// Elastic Arms vertex finding algorithm
namespace earms {
  ///
  /// A stand-alone module for running the elastic arm algorithm
  ///
  class ElasticArmsHS : public art::EDProducer 
  {
  public:
    explicit ElasticArmsHS(fhicl::ParameterSet const &p);
    virtual ~ElasticArmsHS();
    
    void  beginJob();
    void  reconfigure(const fhicl::ParameterSet& p);
    void  produce(art::Event &e);
    art::PtrVector<rb::CellHit> Scrub(art::PtrVector<rb::CellHit> c);
    void  endJob();
    
  private:
    void FindSeed(double xmin,
                  double xmax,
                  double ymin,
                  double ymax,
                  double zmin,
                  double zmax,
                  const rb::HoughResult& hrx,
                  const rb::HoughResult& hry,
                  ElasticArms& arms);
    void MCVertex(art::Event&          evt,
                  std::vector<double>& mcx,
                  std::vector<double>& mcy,
                  std::vector<double>& mcz);
                  
  private:
    int          fCompareToMC;   ///< Should we make a comparison to truth?
    int          fMakeSummaryNt; ///< Make a summary ntuple entry for each event
    int          fMakeHitNt;     ///< Make debugging hit-level ntuple
    unsigned int fMinHit;        ///< Slices must have at least this many hits
    unsigned int fMinHitX;       ///< Miniumm hits in x view / slice
    unsigned int fMinHitY;       ///< Minimum hit in y view / slice
    double       fLambdaV;       ///< Distance scale vertex to first hit
    //
    // Paramters which control the initial grid search
    //
    bool                fBeamBias;           ///< True = assume +z beam bias
                                             ///< False = assume -z beam bias
    double              fZTruncFracLo;       ///< z hit to base the window on (fraction)
    double              fZTruncFracHi;       ///< z hit to base the window on (fraction)
    double              fZTruncDz1;          ///< How far up stream to set window (cm)
    double              fZTruncDz2;          ///< How far down stream to set window (cm)
    double              fVtxBoxBuffer;       ///< Vertex must be this close to hits [cm]
    std::vector<double> fGridSearchFz;       ///< Fractional locations of z vtx seeds
    double              fGridSearchColinear; ///< Colinearity cut
    unsigned int        fGridSearchNdir;     ///< Number of minimum bias directions to use
    double              fGridSearchT;        ///< Temperature to use for grid search
    double              fGridSearchLambda;   ///< Noise penalty for grid search
    unsigned int        fMaxHoughResults;    ///< Only use the top hough results
    unsigned int        fNhoughMatch;        ///< How many lines to check in the other view
    double              fHoughPeakAsymCut;   ///< Which x/y hough results to pair?
    double              fHoughPeakThr;       ///< Threshold to accept a peak as good
    
    
    //
    // Parameters which control the final minimization
    //
    double       fMinimizeLambda;  ///< Noise penalty during minimization
    double       fMinimizeT0;      ///< Initial temperature for annealing
    double       fMinimizeT1;      ///< Final temperature for annealing
    double       fAnnealFrac;      ///< Fractional decrease in temperature
    unsigned int fAnnealMaxItr;    ///< Maximum number of annealing steps
    //
    // Diagnostics
    //
    TNtuple* fSumNt;               ///< Summary ntuple
    TNtuple* fHitNt;               ///< Debugging hit-level ntuple
    TH1F*    fVtxDx;               ///< Distribution of x(fit)-x(mc)
    TH1F*    fVtxDy;               ///< Distribution of y(fit)-y(mc)
    TH1F*    fVtxDz;               ///< Distribution of z(fit)-z(mc)
    
    //folder labels
    const art::ProductToken<std::vector<rb::Cluster>> fSliceToken;
    std::string fHoughResultLabel; ///< Module label to output HoughResults
    std::string fMCLabel;          ///< Module label for MC
    
    bool fObeyPreselection;        ///< Check rb::IsFiltered?
  };
}

namespace earms{
  
  //......................................................................
  
  ElasticArmsHS::ElasticArmsHS(fhicl::ParameterSet const &p) : 
    fHitNt(0),
    fVtxDx(0),
    fVtxDy(0),
    fVtxDz(0),
    fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SliceLabel")))
  {
    this->reconfigure(p);
    this->produces< std::vector<rb::Vertex> >();
    this->produces< std::vector<rb::Prong>  >();
    this->produces<art::Assns<rb::Vertex, rb::Cluster> >();
    this->produces<art::Assns<rb::Prong, rb::Cluster> >();
  }
  
  //......................................................................
  
  void ElasticArmsHS::beginJob() 
  {
    art::ServiceHandle<art::TFileService> f;
    
    if(fMakeSummaryNt) {
      fSumNt = f->make<TNtuple>("nt",
                                "nt",
                                "run:subrun:evt:e:nhit:m:sumvia:xf:yf:zf:xmc:ymc:zmc");
    }
    if(fMakeHitNt) {
      fHitNt = f->make<TNtuple>("hnt","hnt","run:subrun:evt:xy:z:view:a:mia:via");
    }
    if(fCompareToMC) {
      fVtxDx = f->make<TH1F>("fVtxDx",
                             ";x(fit) - x(mc) [cm];Events",
                             100,-50.0,50.0);
      fVtxDy = f->make<TH1F>("fVtxDy",
                             ";y(fit) - y(mc) [cm];Events",
                             100,-50.0,50.0);
      fVtxDz = f->make<TH1F>("fVtxDz",
                             ";z(fit) - z(mc) [cm];Events",
                             100,-50.0,50.0);
    }
  }
  
  //......................................................................
  
  void ElasticArmsHS::reconfigure(const fhicl::ParameterSet& p) 
  {
#define SET(N) { f##N = p.get<__typeof__(f##N)>(#N); }
    
    SET(CompareToMC);
    SET(MakeSummaryNt);
    SET(MakeHitNt);
    
    SET(MinHit);
    SET(MinHitX);
    SET(MinHitY);
    SET(MaxHoughResults);
    SET(NhoughMatch);
    SET(AnnealMaxItr);
    SET(GridSearchNdir);
    
    SET(LambdaV);
    SET(ZTruncFracLo);
    SET(ZTruncFracHi);
    SET(ZTruncDz1);
    SET(ZTruncDz2);
    SET(VtxBoxBuffer);
    SET(HoughPeakAsymCut);
    SET(GridSearchColinear);
    SET(GridSearchT);
    SET(GridSearchLambda);
    SET(MinimizeLambda);
    SET(MinimizeT0);
    SET(MinimizeT1);
    SET(AnnealFrac);
    SET(HoughPeakThr);
    
    SET(HoughResultLabel);
    SET(MCLabel);
    
    SET(ObeyPreselection);
    SET(BeamBias);
    
    SET(GridSearchFz);
#undef SET
  }
  
  //......................................................................
  
  ElasticArmsHS::~ElasticArmsHS() { }
  
  //......................................................................
  
  void ElasticArmsHS::FindSeed(double xmin,
                               double xmax,
                               double ymin,
                               double ymax,
                               double zmin,
                               double zmax,
                               const rb::HoughResult& hrx,
                               const rb::HoughResult& hry,
                               ElasticArms& arms)
  {
    unsigned int i, j;
    GridSearch gs(xmin-fVtxBoxBuffer,
                  xmax+fVtxBoxBuffer,
                  ymin-fVtxBoxBuffer,
                  ymax+fVtxBoxBuffer,
                  zmin-fVtxBoxBuffer,
                  zmax+fVtxBoxBuffer,
                  fGridSearchColinear,
                  fGridSearchNdir);
    //
    // Produce a list of seeds using all intersections of hough lines
    // and points near hits in the event
    //
    std::vector<double> dxv;
    std::vector<double> dyv;
    std::vector<double> bxv;
    std::vector<double> byv;
    for(i = 0; i < fNhoughMatch; ++i) {
      double tmpd, tmpb;
      if(i < hrx.fPeak.size()) {
        hrx.SlopeIcept(i,&tmpd,&tmpb);
        dxv.push_back(tmpd);
        bxv.push_back(tmpb);
      }
      if(i < hry.fPeak.size()) {
        hry.SlopeIcept(i,&tmpd,&tmpb);
        dyv.push_back(tmpd);
        byv.push_back(tmpb);
      }
    }
    gs.AddVtxPoints(arms, dxv, bxv, dyv, byv, fGridSearchFz);
    gs.AddHoughIntersections(hrx, hry, fMaxHoughResults, fNhoughMatch);
    
    //
    // Set the seed directions. Use a standard set plus combinations of
    // the strongest hough lines.
    //
    gs.SetStandardDirections();
    unsigned int n = fMaxHoughResults;
    if(hrx.fPeak.size() < n) n = hrx.fPeak.size();
    if(hry.fPeak.size() < n) n = hry.fPeak.size();
    double hx, hy;
    double dx, bx, dy, by;
    for(i = 0; i < n; ++i) {
      hrx.SlopeIcept(i,&dx,&bx);
      hx = hrx.fPeak[i].fH;
      for(j = 0; j < n; ++j) {
        hry.SlopeIcept(j,&dy,&by);
        hy = hry.fPeak[j].fH;
        //
        // To beat down the combinatorics, only match lines that have
        // nearly the same peak height in each view
        //
        if(fabs(hx-hy)/(hx+hy) < fHoughPeakAsymCut) {
          gs.AddDirection(acos( 1/sqrt(dx*dx+dy*dy+1)),atan2(dy,dx));
          gs.AddDirection(acos(-1/sqrt(dx*dx+dy*dy+1)),atan2(dy,dx));
        }
      }
    }
    
    arms.SetT     (fGridSearchT);
    arms.SetLambda(fGridSearchLambda);
    gs.FindBestVtx(arms);
  }
  
  //....................................................................
  // Find the primary vertex location using true Monte Carlo
  // information. In principle there can be more than one MC
  // interaction, so return them as lists of x,y,z positions.
  //
  void ElasticArmsHS::MCVertex(art::Event&          evt,
                               std::vector<double>& mcx,
                               std::vector<double>& mcy,
                               std::vector<double>& mcz)
                
  {
    unsigned int j;
    art::Handle< std::vector<simb::MCTruth> > mclist;
    evt.getByLabel(fMCLabel, mclist);
    mcx.resize(mclist->size());
    mcy.resize(mclist->size());
    mcz.resize(mclist->size());
    for (j=0; j<mclist->size(); ++j) {
      art::Ptr<simb::MCTruth> mc(mclist,j);
      const simb::MCNeutrino& ne(mc->GetNeutrino());
      const simb::MCParticle& nu(ne.Nu());
      mcx[j] = nu.Vx();
      mcy[j] = nu.Vy();
      mcz[j] = nu.Vz();
    }
  }

  //......................................................................

  void ElasticArmsHS::produce(art::Event& evt) 
  {
    unsigned int i, j;
    
    std::unique_ptr< std::vector<rb::Vertex> > vtx(new std::vector<rb::Vertex>);
    std::unique_ptr< std::vector<rb::Prong>  > prg(new std::vector<rb::Prong>);
    std::unique_ptr< art::Assns<rb::Vertex, rb::Cluster> > 
      vAssns(new art::Assns<rb::Vertex, rb::Cluster>);
    std::unique_ptr< art::Assns<rb::Prong, rb::Cluster> > 
      pAssns(new art::Assns<rb::Prong, rb::Cluster>);
    
    art::ServiceHandle<geo::Geometry> fGeo;
    
    //
    // Pull out the time slices and hough results for the event
    //
    art::Handle<std::vector<rb::Cluster> > slice;
    evt.getByToken(fSliceToken,slice);

    // Object to get hough lines associated with a slice
    art::FindMany<rb::HoughResult> fmhr(slice, evt, fHoughResultLabel);
    
    for(i = 0; i < slice->size(); ++i) {
      if(fObeyPreselection && rb::IsFiltered(evt, slice, i)) continue;
      const rb::Cluster& s = (*slice)[i];
      if(s.IsNoise()) continue;
      // Scrub this cluster of isolated hits.
      // When slicer is fixed this will no longer be necessary
      rb::Cluster neighbors(this->Scrub(s.AllCells()));
      //
      // Skip the noise cluster and require some minimum number of hits
      //
      if(neighbors.NXCell() < fMinHitX) continue;
      if(neighbors.NYCell() < fMinHitY) continue;
      if(neighbors.NCell()  < fMinHit)  continue;
      
      //
      // Find any hough results that go with this time slice
      //
      const rb::HoughResult* hrx = 0;
      const rb::HoughResult* hry = 0;
      
      std::vector<const rb::HoughResult*>  hr = fmhr.at(i);
      for(j = 0; j < hr.size(); ++j) {
        if(hr[j]->fView == geo::kX) {
          hrx = hr[j];
        }
        if(hr[j]->fView == geo::kY) {
          hry = hr[j];
        }
      }
      //
      // No hough results, no joy.
      //
      if(hrx == 0 || hry == 0)  continue;
      if(hrx->fPeak.size() < 1) continue;
      if(hry->fPeak.size() < 1) continue;
      
      //
      // Decide how many prongs we have
      //
      unsigned int nx = 0;
      for(j = 0; j < hrx->fPeak.size(); ++j) {
        if(hrx->fPeak[j].fH/hrx->fPeak[j].fThr >= fHoughPeakThr) ++nx;
      }
      unsigned int ny = 0;
      for(j = 0; j < hry->fPeak.size(); ++j) {
        if(hry->fPeak[j].fH/hry->fPeak[j].fThr >= fHoughPeakThr) ++ny;
      }
      unsigned int nprong = 1;
      if(nx > nprong) nprong = nx;
      if(ny > nprong) nprong = ny;
      if(nprong > fMaxHoughResults) nprong = fMaxHoughResults;
      
      //
      // Decide where to truncate the event along z
      //
      std::vector<double> zs;
      for (j = 0; j < neighbors.NCell(); ++j) {
        const rb::CellHit* h = neighbors.Cell(j).get();
        
        unsigned short p = h->Plane();
        unsigned short c = h->Cell();
        
        double      xyz[3];
        double      dxyz[3];
        geo::View_t view;
        fGeo->CellInfo(p, c, &view, xyz, dxyz);
        
        // To form a vertex position that has all the implicit upstream biased
        // reversed, just flip all the points around before examining them.
        zs.push_back(fBeamBias ? xyz[2] : -xyz[2]);
      }
      sort(zs.begin(), zs.end());

      unsigned int iz = (fZTruncFracLo*zs.size())/100;
      double zlo = zs[iz] - fZTruncDz1;
      double zhi = zs[iz] + fZTruncDz2;

      unsigned int ncell = 0, ic = 0;
      std::vector<double> zstemp;
      for(j = 0; j < zs.size(); ++j) {
        if(zs[j] >= zlo && zs[j] <= zhi){
          ++ncell;
          zstemp.push_back(zs[j]);
        }
      }

      if(ncell == 0) { 
        zstemp.clear();
        zlo = zs[0]           - fZTruncDz1; 
        zhi = zs[zs.size()-1] + fZTruncDz2;
        for(j = 0; j < zs.size(); ++j) {
          if(zs[j] >= zlo && zs[j] <= zhi){
            ++ncell;
            zstemp.push_back(zs[j]);
          }
        }
      }
      
      //
      // Construct the arms and fill with cell hits
      //
      unsigned int nhitx = 0, nhity = 0;
      ElasticArms arms(ncell, nprong, fGridSearchLambda, fLambdaV);
      double minZ = fBeamBias ? fGeo->DetLength() : 0;
      double maxZ = fBeamBias ? 0 : fGeo->DetLength();
      double minX = fGeo->DetHalfWidth();
      double maxX = -fGeo->DetHalfWidth();
      double minY = fGeo->DetHalfHeight();
      double maxY = -fGeo->DetHalfHeight();
      for(j = 0; j < neighbors.NCell(); ++j) {
        const rb::CellHit* h = neighbors.Cell(j).get();
        
        unsigned short p = h->Plane();
        unsigned short c = h->Cell();
        
        double      xyz[3];
        double      dxyz[3];
        geo::View_t view;
        fGeo->CellInfo(p, c, &view, xyz, dxyz);

        if(xyz[2] >= zlo && xyz[2] <= zhi) {
          if(view == geo::kX) {
            arms.SetHit(ic++, xyz[2], xyz[0], dxyz[2]/sqrt(12.), view);
            ++nhitx;
            // For the purpose of finding the Hough vertex seed, don't flip
            // anything (ie we have to put it back). This prevents having to edit
            // any of that code.
            if(xyz[2] > maxZ) maxZ = fBeamBias ? xyz[2] : -xyz[2];
            if(xyz[2] < minZ) minZ = fBeamBias ? xyz[2] : -xyz[2];
            if(xyz[0] > maxX) maxX = xyz[0];
            if(xyz[0] < minX) minX = xyz[0];
          }
          else if(view == geo::kY) {
            arms.SetHit(ic++, xyz[2], xyz[1], dxyz[2]/sqrt(12.), view);
            ++nhity;
            if(xyz[2] > maxZ) maxZ = fBeamBias ? xyz[2] : -xyz[2];
            if(xyz[2] < minZ) minZ = fBeamBias ? xyz[2] : -xyz[2];
            if(xyz[1] > maxY) maxY = xyz[1];
            if(xyz[1] < minY) minY = xyz[1];
          }
          else abort();
        }
      }
      
      //
      // Need at least one hit in each view to continue
      //
      if(nhitx == 0 || nhity == 0) continue;
      
      //
      // Seed the arms before trying to minimize
      //
      this->FindSeed(minX, maxX, minY, maxY, minZ, maxZ, *hrx, *hry, arms);

      if(!fBeamBias){
        // Flip the resulting seed around to match the flipped points we'll be
        // using.
        double vx, vy, vz;
        arms.GetVertex(vx, vy, vz);
        arms.SetVertex(vx, vy, -vz);
      }
      
      //
      // With the arms seeded, now try to find the best fit
      //   
      //arms.SetVertex(415.0,-603.0,4980.0);
      arms.SetLambda(fMinimizeLambda);
      Minimizer mn(&arms);
      mn.SetT0(fMinimizeT0);
      mn.SetT1(fMinimizeT1);
      mn.SetAnnealFrac(fAnnealFrac);
      mn.SetAnnealMaxItr(fAnnealMaxItr);
      mn.Fit();
      
      //
      // Fill the hit-level debugging ntuple
      //
      if(fMakeHitNt) {
        float x[9];
        for(j = 0; j < arms.fN; ++j) {
          x[0] = evt.run();
          x[1] = evt.subRun();
          x[2] = evt.id().event();
          x[3] = arms.fXorY[j];
          x[4] = arms.fZ[j];
          x[5] = arms.fView[j];
          for(unsigned int a = 0; a < arms.fM; ++a) {
            x[6] = a;
            x[7] = arms.fMia[j][a];
            x[8] = arms.fVia[j][a];
            fHitNt->Fill(x);
          }
        }
      }
      
      //
      // Push the results out
      //
      art::Ptr<rb::Cluster> sptr(slice, i);
      rb::Vertex v(arms.fX0, arms.fY0, fBeamBias ? arms.fZ0 : -arms.fZ0,
                   neighbors.MinTNS());
      vtx->push_back(v);
      util::CreateAssn(*this, evt, *vtx, sptr, *vAssns);
      
      for(j = 0; j < arms.fdXds.size(); ++j) {
        rb::Prong p(s,
                    TVector3(arms.fX0,      arms.fY0,      arms.fZ0),
                    TVector3(arms.fdXds[j], arms.fdYds[j], arms.fdZds[j]));
        prg->push_back(p);
        util::CreateAssn(*this, evt, *prg, sptr, *pAssns);
      }
      
      //
      // If we've been asked, make a comparison to MC truth
      //
      std::vector<double> xmc; // MC x vertex locations
      std::vector<double> ymc; // MC y vertex locations
      std::vector<double> zmc; // MC z vertex locations
      if (fCompareToMC) {
        this->MCVertex(evt, xmc, ymc, zmc);
        for(j = 0; j < xmc.size(); ++j) {
          double dx = arms.fX0 - xmc[j];
          double dy = arms.fY0 - ymc[j];
          double dz = arms.fZ0 - zmc[j];
          fVtxDx->Fill(dx);
          fVtxDy->Fill(dy);
          fVtxDz->Fill(dz);
          if(fabs(dx) > 4*4 || fabs(dy) > 4*4 || fabs(dz) > 4*6) {
            std::cout << "ElasticArmsHS: Blown fit [" 
                      << evt.run() << ":" 
                      << evt.subRun() << ":" 
                      << evt.id().event() << "] dx,dy,dz=("
                      << dx << ","
                      << dy << ","
                      << dz << ")" 
                      << std::endl;
          } // if we have a blown fit
        } // for each entry in mclist
      } // if making mc comparisons
      
      //
      // Record a summary entry to the ntuple
      //
      if(fMakeSummaryNt) {
        double sum = 0.0;
        for(i = 0; i < arms.fN; ++i) {
          for(j = 0; j < arms.fM; ++j) {
            sum += arms.fVia[i][j];
          }
        }
        float x[13];
        x[ 0] = evt.run();
        x[ 1] = evt.subRun();
        x[ 2] = evt.id().event();
        x[ 3] = arms.E();
        x[ 4] = arms.fN;
        x[ 5] = arms.fM;
        x[ 6] = sum;
        x[ 7] = arms.fX0;
        x[ 8] = arms.fY0;
        x[ 9] = arms.fZ0;
        if (xmc.size()>0) {
          x[10] = xmc[0];
          x[11] = ymc[0];
          x[12] = zmc[0];
        }
        else {
          x[10] = -9e9;
          x[11] = -9e9;
          x[12] = -9e9;
        }
        fSumNt->Fill(x);
      } // if we are making a summary ntuple
  
    } // Overall loop on slices

    evt.put(std::move(vtx));
    evt.put(std::move(prg));
    evt.put(std::move(vAssns));
    evt.put(std::move(pAssns));
  }
  
  //......................................................................
  //function to remove hits with no neighbors from the slice
  art::PtrVector<rb::CellHit> ElasticArmsHS::Scrub(art::PtrVector<rb::CellHit> c)
  {
    
    std::vector<int> hit_numNeighbors;
    //in an ideal detector a hit must have 1 neighbor
    //inside this distance.
    double maxGap = 60.0; 
    art::PtrVector<rb::CellHit> trimCells;
    int iPlane, fPlane, iCell, fCell;
    double goodCh, totCh;
    
    iPlane = 0;
    iCell  = 0;
    fPlane = 0;
    fCell  = 0;
    
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<chaninfo::BadChanList> badc;
    
    for(unsigned int i = 0; i < c.size(); i++){
      hit_numNeighbors.push_back(0);
    }
    
    //loop to find a cells closest neighbor
    for(unsigned int i = 0; i < c.size(); i++){
      double minDist = 9999.0;
      int minCell = -1;
      if(hit_numNeighbors[i] > 0) continue;
      for(unsigned int j = 0; j < c.size(); j++){
        if(i == j) continue;
        if(c[i]->View() != c[j]->View()) continue;
        double xyz1[3], xyz2[3];
        geom->Plane((c[i])->Plane())->Cell((c[i])->Cell())->GetCenter(xyz1);
        geom->Plane((c[j])->Plane())->Cell((c[j])->Cell())->GetCenter(xyz2);
        double dist = util::pythag(xyz1[0]-xyz2[0],
                                   xyz1[1]-xyz2[1],
                                   xyz1[2]-xyz2[2]);
        if(dist < minDist){
          minDist = dist;
          minCell = j;
        }
      }
      if(minCell == -1) continue;
      
      //define the rectangle of cells such that a cell and its 
      //nearest neighbor are at diagonal corners
      if(c[i]->Plane() < c[minCell]->Plane()){
        iPlane = c[i]->Plane();
        fPlane = c[minCell]->Plane();
      }
      else {
         iPlane = c[minCell]->Plane();
         fPlane = c[i]->Plane();
      }
      if(c[i]->Cell() < c[minCell]->Cell()) {
        iCell = c[i]->Cell();
        fCell = c[minCell]->Cell();
      }
      else {
        iCell = c[minCell]->Cell();
        fCell = c[i]->Cell();
      }
      
      //within the rectangle find the fraction of cells that are good
      goodCh = 0;
      totCh = 0;
      for(int c = iCell; c <= fCell; c++){
        for(int p = iPlane; p <= fPlane; p+=2){
          totCh++;
          if(!(badc->IsBad(p,c))) ++goodCh;
        }
      }
      if((minDist*goodCh/totCh) < maxGap){
        hit_numNeighbors[i]++;
        hit_numNeighbors[minCell]++;
      }
    }
    
    //keep the hits that have at least 1 neighbor in the slice
    for(unsigned int i = 0; i < c.size(); i++){
      if(hit_numNeighbors[i] > 0){
        trimCells.push_back(c[i]);
      }
    }
   
    return trimCells;
  }
  
  
  //......................................................................
  
  void ElasticArmsHS::endJob() 
  {
    if(fCompareToMC) {
      std::cout << "ElasticArmsHS: dx | rms=" << fVtxDx->GetRMS() 
                << " | outliers=" 
                << fVtxDx->Integral(0,30)+fVtxDx->Integral(70,101) 
                << std::endl;
      std::cout << "ElasticArmsHS: dy | rms=" << fVtxDy->GetRMS() 
                << " | outliers=" 
                << fVtxDy->Integral(0,30)+fVtxDy->Integral(70,101) 
                << std::endl;
      std::cout << "ElasticArmsHS: dz | rms=" << fVtxDz->GetRMS() 
                << " | outliers=" 
                << fVtxDz->Integral(0,30)+fVtxDz->Integral(70,101) 
                << std::endl;
    }
  }

}
  ////////////////////////////////////////////////////////////////////////
namespace earms
{
  DEFINE_ART_MODULE(ElasticArmsHS)
}
////////////////////////////////////////////////////////////////////////