BreakPoint_module.cc

Go to the documentation of this file.

///
/// \file    BreakPoint_module.cc
/// \brief   Run the break-point fitter
/// \author  messier@indiana.edu
/// \version $Id:$
///
#include "TH2F.h"
#include "TNtuple.h"
//
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
//
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Track.h"
#include "RecoBase/FitSum.h"
//
#include "Utilities/AssociationUtil.h"
//
#include "BreakPointFitter/func/HitList3D.h"
#include "BreakPointFitter/func/TrackBasis.h"
#include "BreakPointFitter/func/Path.h"
#include "BreakPointFitter/func/ScatteringSurfaces.h"
#include "BreakPointFitter/func/Lutz.h"
#include "BreakPointFitter/func/BPException.h"

namespace bpfit { class BreakPoint; }

class bpfit::BreakPoint : public art::EDProducer
{
public:
  const art::ProductToken<std::vector<rb::Cluster>> fSliceToken;

  // FHICL Configuration data
  bool         fDebugHitList0;   ///< Debug the detector-frame hit list?
  bool         fDebugHitList1;   ///< Debug the track-frame hit list?
  bool         fDebugScatSurf;   ///< Debug the scattering surfaces
  bool         fDebugFitList0;   ///< Debug the fitted track locations (det. frame)
  bool         fDebugFitList1;   ///< Debug the fitted track locations (track frame)
  bool         fLogResiduals;    ///< Make ntuple of residuals?
  bool         fFillTrackNt;     ///< Fill ntuple of tracking results?
  bool         fFillExceptNt;    ///< Fill ntuple of exception cases?
  bool         fVerbose;         ///< Run in verbose mode (print messages about all exceptions?)
  bool         fFitAsMuon;       ///< Fit tracks under the assumption of a muon?
  bool         fFitAsPion;       ///< Fit tracks under the assumption of a pion?
  bool         fFitAsProton;     ///< Fit tracks under the assumption of a proton?
  double       fExtendUp;        ///< Fractionally, how far to extend track into upstream dead material
  double       fExtendDn;        ///< Fractionally, how far to extend track into downstream dead material
  std::string  fVertexLabel;     ///< Where to find vertices in event
  std::string  fClusterLabel;    ///< Where to find clusters in event
  std::string  fClusterInstance; ///< Which instance of clusters to find
  unsigned int fNhitX;           ///< Require prongs to have this many x hits
  unsigned int fNhitY;           ///< Require prongs to have this many y hits
  unsigned int fNhit;            ///< Require this many hits total
  double       fChi2Cutoff;      ///< Drop hits with this chi2 or greater
  double       fdStol;           ///< Scattering surface ds tolerance
  double       fdXtol;           ///< Scattering surface dx tolerance
  double       fdTtol;           ///< Scattering surface dt tolerance
  double       fSigAdj;          ///< Multiplier to adjust sigmaXYZ for hits (see comment in the FillHitList function.) Set to 1.0 if no adjustment is desired.
  unsigned int fMaxSliceHits;    ///< Don't bother processing a slice if it has more than this many hits. This is to prevent computing problems.

public:
  // Member data
  int          fRun;
  int          fSubrun;
  int          fEvent;
  unsigned int fSlice;
  unsigned int fProng;
  TNtuple*     fHitList0Nt;
  TNtuple*     fHitList1Nt;
  TNtuple*     fScatSurfNt;
  TNtuple*     fFitList0Nt;
  TNtuple*     fFitList1Nt;
  TNtuple*     fResidNt;
  TNtuple*     fTrackNt;
  TNtuple*     fExceptNt;

  std::vector<int> fPartList;    ///< List of particle types (by pdg code) to be used for track fits.

public:
  //......................................................................
  
  void reconfigure(fhicl::ParameterSet const& p)
  {
    fDebugHitList0 = p.get<bool>("DebugHitList0");
    fDebugHitList1 = p.get<bool>("DebugHitList1");
    fDebugScatSurf = p.get<bool>("DebugScatSurf");
    fDebugFitList0 = p.get<bool>("DebugFitList0");
    fDebugFitList1 = p.get<bool>("DebugFitList1");
    fLogResiduals  = p.get<bool>("LogResiduals");
    fFillTrackNt   = p.get<bool>("FillTrackNt");
    fFillExceptNt  = p.get<bool>("FillExceptNt");
    fVerbose       = p.get<bool>("Verbose");

    fExtendUp    = p.get<double>("ExtendUp");
    fExtendDn    = p.get<double>("ExtendDn");
    fFitAsMuon   = p.get<bool>("FitAsMuon");
    fFitAsPion   = p.get<bool>("FitAsPion");
    fFitAsProton = p.get<bool>("FitAsProton");
    fPartList.clear();
    if(fFitAsMuon)   fPartList.push_back(13);
    if(fFitAsPion)   fPartList.push_back(211);
    if(fFitAsProton) fPartList.push_back(2212);

    fVertexLabel     = p.get<std::string>("VertexLabel");
    fClusterLabel    = p.get<std::string>("ClusterLabel");
    fClusterInstance = p.get<std::string>("ClusterInstance");

    fNhitX      = p.get<unsigned int>("NhitX");
    fNhitY      = p.get<unsigned int>("NhitY");
    fNhit       = p.get<unsigned int>("Nhit");
    fChi2Cutoff = p.get<double>("Chi2Cutoff");
    fdStol      = p.get<double>("dStol");
    fdXtol      = p.get<double>("dXtol");
    fdTtol      = p.get<double>("dTtol");
    fSigAdj     = p.get<double>("SigmaAdjust");
    fMaxSliceHits = p.get<unsigned int>("MaxSliceHits");
  }
  
  //......................................................................

  explicit BreakPoint(fhicl::ParameterSet const& p) :
    fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SliceLabel"))),
    fRun(0),
    fSubrun(0),
    fEvent(0),
    fProng(0),
    fHitList0Nt(0),
    fHitList1Nt(0),
    fScatSurfNt(0),
    fFitList0Nt(0),
    fFitList1Nt(0),
    fResidNt(0),
    fTrackNt(0),
    fExceptNt(0)
  {
    produces< std::vector<rb::Track> >();
    produces< std::vector<rb::FitSum> >();
    produces< art::Assns<rb::Track, rb::Prong> >();
    produces< art::Assns<rb::Track, rb::Cluster> >();
    produces< art::Assns<rb::FitSum, rb::Track> >();
    this->reconfigure(p);
  }

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

  virtual ~BreakPoint() { }

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

  void FindSlices(const art::Event& evt, 
                  art::Handle< std::vector<rb::Cluster> >& h,
                  art::PtrVector<rb::Cluster>& s) 
  {
    unsigned int i;
    evt.getByToken(fSliceToken, h);
    for (i=0; i<h->size(); ++i) {
      s.push_back(art::Ptr<rb::Cluster>(h, i));
    }
  }

  //......................................................................
  
  void FindProngs(art::Event&                          evt,
                  unsigned int                         islice,
                  art::FindManyP<rb::Prong>&           prongh,
                  std::vector< art::Ptr<rb::Vertex> >& vertexp,
                  std::vector< art::Ptr<rb::Prong> >&  prongp)
  {
    if (prongh.isValid()) {
      prongp = prongh.at(islice);
    }
    else {
      //
      // Hmm... No prongs associated with slice - try the vertex.
      //
      art::FindManyP<rb::Prong> prongh2(vertexp, evt, fClusterLabel);
      if (prongh2.isValid()) {
        if (prongh2.size()==1) {
          prongp = prongh2.at(0);
        }
        else {
          std::cerr << __FILE__ << ":" << __LINE__ 
                    << " Multiple vertices breaks our assumptions."
                    << std::endl;
        }
      }
    }
  }
  
  //......................................................................
  
  bool GoodProng(art::Ptr<rb::Prong>& t)
  {
    if (t->Is2D()) return false;
    
    unsigned int nx = t->NXCell();
    unsigned int ny = t->NYCell();
    if (nx<fNhitX || ny<fNhitY) return false;
    if (nx+ny<fNhit)            return false;
    
    return true;
  }
  
  //......................................................................
  
  void FillHitList(art::Ptr<rb::Prong>& prongp,
                   bpfit::HitList3D&    hits)
  {
    art::ServiceHandle<geo::Geometry> geom;
    unsigned int itr, i;
    unsigned int n = prongp->NCell();
    //
    // Build the 3D hits in two passes. In the first we don't have
    // information about where transversly we are along a cell. In the
    // second iteration, we do.
    //
    for (itr=0; itr<2; ++itr) {
      for (i=0; i<n; ++i) {
        const rb::CellHit* h = prongp->Cell(i).get();
        
        unsigned int ip = h->Plane();
        unsigned int ic = h->Cell();
        
        const geo::PlaneGeo* p = geom->Plane(ip);
        const geo::CellGeo*  c = p->   Cell(ic);
        
        double xyz[3] = {0,0,0};;
        double sigx=0, sigy=0, sigz=0, sigt=50;
        double w = 1.0;
        if (itr==0) {
          //
          // On first iteration we don't have information from other
          // view. Use cell centers for position estimates.
          //

          // NOTE: BPF appears to be giving too much weight to the hits
          //       when doing the track fit (the total chi2 distribution
          //       does not have a mean of 1.0.) I believe this to be
          //       related to the fact that the true particle location
          //       within a cell is not Gaussian distributed around the
          //       cell center. Introducing the fSigAdj factor adjusts
          //       this and brings the mean of the chi2 distribution down
          //       to 1.0 allowing for the fit scattering angles to be better.
          c->GetCenter(xyz);
          if (p->View()==geo::kX) sigx = c->HalfW()*fSigAdj;
          if (p->View()==geo::kY) sigy = c->HalfW()*fSigAdj;
          sigz  = c->HalfD()*fSigAdj;
        }
        else {
          //
          // On second iteration, we have more information about where
          // in the cell we are. Use the estimated location from the
          // other view. Otherwise, most information is a copy of what
          // is already there.
          //
          if (p->View()==geo::kX) c->GetCenter(xyz, hits[i].fY);
          if (p->View()==geo::kY) c->GetCenter(xyz, hits[i].fX);
          sigx = hits[i].fSigX;
          sigy = hits[i].fSigY;
          sigz = hits[i].fSigZ;
          sigt = hits[i].fSigT;
          w    = hits[i].fW;
        }
        
        if (p->View()==geo::kX) {
          hits.SetHit(i, geo::kX, xyz[0], xyz[1], xyz[2], h->TNS(),
                      sigx, sigy, sigz, sigt, w);
        }
        else if (p->View()==geo::kY) {
          hits.SetHit(i,
                      geo::kY, xyz[0], xyz[1], xyz[2], h->TNS(),
                      sigx, sigy, sigz, sigt, w);
        }
        else abort();
      } // loop on hits
      hits.CompleteOrtho();
    } // iterations
    //
    // Keep all hits inside the detector
    //
    double xyz[3];
    for (i=0; i<hits.size(); ++i) {
      xyz[0] = hits[i].fX;
      xyz[1] = hits[i].fY;
      xyz[2] = hits[i].fZ;
      this->ClampToDetector(xyz);
      hits[i].fX = xyz[0];
      hits[i].fY = xyz[1];
      hits[i].fZ = xyz[2];
    }
    
    hits.SortByZ();
  }

  //......................................................................
  ///
  /// If there are less than 10 hits in one view with the same z then 
  /// take the average x or y depending on view and use that when 
  /// fitting the track.
  ///
  void AverageXY(bpfit::HitList3D& h)
  {

    // NOTE: To the reader of this code - I think there is a mistake below.
    //       I don't think the intention is to set the weight of any hits
    //       to 1e-9. I think that the weights should be left as 1.0. But
    //       we aren't currently using this function...

    for(unsigned int i=0;i<h.size();i++){
      
      int lastj = i;
      int n_xhits = 0;
      int n_yhits = 0;
      unsigned int i_xhits[384], i_yhits[384]; // Max number of hits at the same z = # cells per plane
      double  sumx = 0.0;
      double  sumy = 0.0;
      double sigx2 = 0.0;
      double sigy2 = 0.0;
      
      for(unsigned int j=i+1;j<h.size();j++){   
        // Look for hits in the same view that have the 
        // same z coordinate.
        if((fabs(h[i].fZ-h[j].fZ)<0.1)){
          if(h[i].fView == h[j].fView){
            lastj = j;
            if(h[i].fView == geo::kX){   
              if(n_xhits==0) i_xhits[n_xhits] = i;
              else i_xhits[n_xhits] = j;
              n_xhits++;
            }
            else{
              if(h[i].fView == geo::kY){
                if(n_yhits==0) i_yhits[n_yhits] = i;
                else i_yhits[n_yhits] = j;
                n_yhits++;
              }
            }
          }
        }
        else continue;  
      } // end j loop over hitlist3D

      if(n_xhits>0 && n_xhits<10){
        // Averaging x
        for(int xi=0; xi<n_xhits; xi++){
          sumx  += h[i_xhits[xi]].fX;
          sigx2 += h[i_xhits[xi]].fSigX*h[i_xhits[xi]].fSigX;

          // Since we will be adding the averaged hit to 
          // the hitlist we need to reweight the same z hits
          // to something very small.
          if(xi!=0)h[i_xhits[xi]].fW = 1.0E-9;
        }
        double xave = sumx/n_xhits;
        // Replace first hit in series that has the same z with 
        // new average x coordinate.
        h.SetHit(i_xhits[0], h[i_xhits[0]].fView, xave, h[i_xhits[0]].fY, h[i_xhits[0]].fZ, 
                 h[i_xhits[0]].fT, sqrt(sigx2), h[i_xhits[0]].fSigY,
                 h[i_xhits[0]].fSigZ, h[i_xhits[0]].fSigT, h[i_xhits[0]].fW);   
      } // end n_xhits requirement

      if(n_yhits>0 && n_yhits<10){
        // Averaging y
        for(int yi=0; yi<n_yhits; yi++){
          sumy  += h[i_yhits[yi]].fY;
          sigy2 += h[i_yhits[yi]].fSigY*h[i_yhits[yi]].fSigY;
          // Since we will be adding the averaged hit to 
          // the hitlist we need to reweight the same z hits
          // to something very small.
          if(yi!=0) h[i_yhits[yi]].fW = 1.0E-9;
        }
        double yave = sumy/n_yhits;   
        // Replace first hit in series that has the same z with 
        // new average y coordinate.
        h.SetHit(i_yhits[0], h[i_yhits[0]].fView, h[i_yhits[0]].fX, yave, h[i_yhits[0]].fZ,
                 h[i_yhits[0]].fT, h[i_yhits[0]].fSigX, sqrt(sigy2),
                 h[i_yhits[0]].fSigZ, h[i_yhits[0]].fSigT, h[i_yhits[0]].fW);
      } // end n_yhits requirement 

      i = lastj;
    } // end i loop over hitlist3D

    h.CompleteOrtho();

    h.SortByZ();

  } // end function definition 
            
  //......................................................................
  ///
  /// Dump the hit list to an ntuple we can study offline for
  /// debugging
  ///
  void DebugHitList(bpfit::HitList3D& h, int which) 
  {
    TNtuple* nt = 0;

    if (which==0) {
      if (fHitList0Nt==0) {
        art::ServiceHandle<art::TFileService> tfs;
        fHitList0Nt = tfs->make<TNtuple>("hits0",
                                         "Det Frame Hit List Debugging",
                                         "run:subrun:evt:trk:v:x:y:z:sx:sy:sz:slice");
      }
      nt = fHitList0Nt;
    }
    if (which==1) {
      if (fHitList1Nt==0) {
        art::ServiceHandle<art::TFileService> tfs;
        fHitList1Nt = tfs->make<TNtuple>("hits1",
                                         "Trk Frame Hit List Debugging",
                                         "run:subrun:evt:trk:v:x:y:z:sx:sy:sz:slice");
      }
      nt = fHitList1Nt;
    }
    float x[12];
    x[0] = fRun;
    x[1] = fSubrun;
    x[2] = fEvent;
    x[3] = fProng;
    x[11]= fSlice;
    for (unsigned int i=0; i<h.size(); ++i) {
      x[4]  = h[i].fView;
      x[5]  = h[i].fX;
      x[6]  = h[i].fY;
      x[7]  = h[i].fZ;
      x[8]  = h[i].fSigX;
      x[9]  = h[i].fSigY;
      x[10] = h[i].fSigZ;
      nt->Fill(x);
    }
  }

  //......................................................................
  
  void FillTrackNt(const HitList3D& h,
                   const Lutz& lx,
                   const Lutz& ly,
                   const Path& path,
                   int pdg_id)
  {
    if (fTrackNt==0) {
      art::ServiceHandle<art::TFileService> tfs;
      fTrackNt = tfs->make<TNtuple>("trks",
                                    "Track summary",
                                    "run:subrun:evt:slice:trk:"
                                    "nhitx:nhity:x0:y0:z0:x1:y1:z1:length:p:"
                                    "theta:phi:"
                                    "chi2xh:chi2xa:chi2yh:chi2ya:pdgID");
    }
    unsigned int i;
    unsigned int nhitx = 0, nhity = 0;
    for (i=0; i<h.size(); ++i) {
      if (h[i].fView==geo::kX && h[i].fW > 1.0E-9) ++nhitx;
      if (h[i].fView==geo::kY && h[i].fW > 1.0E-9) ++nhity;
    }
    
    double p, len;
    len = path.Length();
    if      (pdg_id == 13)   path.MuonParams(0, 0, &p, 0, 0);
    else if (pdg_id == 211)  path.PionParams(0, 0, &p, 0, 0);
    else if (pdg_id == 2212) path.ProtonParams(0, 0, &p, 0, 0);
    // Can't get here unless we are one of the above particle types,
    // but why leave traps...?
    else                     path.MuonParams(0, 0, &p, 0, 0);
    i = h.size()-1;
    double dx = path.fXYZ[1].X() - path.fXYZ[0].X();
    double dy = path.fXYZ[1].Y() - path.fXYZ[0].Y();
    double dz = path.fXYZ[1].Z() - path.fXYZ[0].Z();
    if (dz==0.0) dz = 1.0E-9;
    
    float x[22];
    x[0]  = fRun;
    x[1]  = fSubrun;
    x[2]  = fEvent;
    x[3]  = fSlice;
    x[4]  = fProng;
    x[5]  = nhitx;
    x[6]  = nhity;
    x[7]  = h[0].fX;
    x[8]  = h[0].fY;
    x[9]  = h[0].fZ;
    x[10] = h[i].fX;
    x[11] = h[i].fY;
    x[12] = h[i].fZ;
    x[13] = len;
    x[14] = p;
    x[15] = 180.0*acos(dz/sqrt(dx*dx+dy*dy+dz*dz))/M_PI;
    x[16] = 180.0*atan2(dy,dx)/M_PI;
    x[17] = lx.Chi2XI();
    x[18] = lx.Chi2Beta();
    x[19] = ly.Chi2XI();
    x[20] = ly.Chi2Beta();
    x[21] = pdg_id;
    fTrackNt->Fill(x);
  }

  //......................................................................
  
  void FillExceptNt(BPException Ex)
  {

   if (fExceptNt == 0) {
      art::ServiceHandle<art::TFileService> tfs;
      fExceptNt = tfs->make<TNtuple>("except",
                                    "Exception summary",
                                    "run:subrun:evt:slice:prong:flag:value");
   }

   float x[7];
   x[0]  = fRun;
   x[1]  = fSubrun;
   x[2]  = fEvent;
   x[3]  = fSlice;
   x[4]  = fProng;
   x[5]  = Ex.fFlag;
   x[6]  = Ex.fValue;
   fExceptNt->Fill(x);

  }

  //......................................................................
  
  void MakeBasis(art::Ptr<rb::Vertex>& vtx,
                 const HitList3D&      hits,
                 bpfit::TrackBasis&    sys) 
  {
    //
    // Extract hits into format expected by TrackBasis
    //
    unsigned int i;
    std::vector<double> x(hits.size());
    std::vector<double> y(hits.size());
    std::vector<double> z(hits.size());
    std::vector<double> dx(hits.size());
    std::vector<double> dy(hits.size());
    std::vector<double> dz(hits.size());

    for (i=0; i<hits.size(); ++i) {
      x[i]  = hits[i].fX;
      y[i]  = hits[i].fY;
      z[i]  = hits[i].fZ;
      dx[i] = hits[i].fSigX;
      dy[i] = hits[i].fSigY;
      dz[i] = hits[i].fSigZ;
    }

    sys.MakeBasis(x,y,z,dx,dy,dz,vtx->GetX(),vtx->GetY(),vtx->GetZ());
  }

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

  void ClampToDetector(double* p)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if (p[0]>geom->DetHalfWidth())  p[0] = geom->DetHalfWidth();
    if (p[1]>geom->DetHalfHeight()) p[1] = geom->DetHalfHeight();
    if (p[2]>geom->DetLength())     p[2] = geom->DetLength();

    if (p[0]<-geom->DetHalfWidth())  p[0] = -geom->DetHalfWidth();
    if (p[1]<-geom->DetHalfHeight()) p[1] = -geom->DetHalfHeight();
    if (p[2]<0)                      p[2] = 0;
  }

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

  void FindEndPoints(HitList3D&         h,
                     bpfit::TrackBasis& sys,
                     double*            p0,
                     double*            p1) 
  {
    unsigned int i;
    double zmin =  9e9;
    double zmax = -9e9;
    for (i=0; i<h.size(); ++i) {
      if (h[i].fZ<zmin) zmin = h[i].fZ;
      if (h[i].fZ>zmax) zmax = h[i].fZ;
    }
    p0[0] = 0.0; p0[1] = 0.0; p0[2] = zmin;
    p1[0] = 0.0; p1[1] = 0.0; p1[2] = zmax;
    double dx;
    double dy;
    double dz;
    sys.TrkToDet(0,      0,      zmin,   0,   0,   0,
                 &p0[0], &p0[1], &p0[2], &dx, &dy, &dz);
    sys.TrkToDet(0,      0,      zmax,   0,   0,   0,
                 &p1[0], &p1[1], &p1[2], &dx, &dy, &dz);
    this->ClampToDetector(p0);
    this->ClampToDetector(p1);
  }
  
  //......................................................................

  void DebugScatSurf(bpfit::ScatteringSurfaces& s, int pdgid)
  {
    if (fScatSurfNt==0) {
      art::ServiceHandle<art::TFileService> tfs;
      fScatSurfNt = tfs->make<TNtuple>("scat",
                                       "Scattering Surfaces",
                                       "run:subrun:evt:trk:s:sigalpha:pdgid:slice");
    }
    float x[8];
    x[0] = fRun;
    x[1] = fSubrun;
    x[2] = fEvent;
    x[3] = fProng;
    x[7] = fSlice;
    for (unsigned int i=0; i<s.N(); ++i) {
      x[4] = s.S(i);
      x[5] = s.SigAlpha(i);
      x[6] = pdgid;
      fScatSurfNt->Fill(x);
    }
  }

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

  void DebugFitList(unsigned int which,
                    const std::vector<double>& x,
                    const std::vector<double>& y,
                    const std::vector<double>& z) 
  {
    TNtuple* nt = 0;
    if (which==0) {
      if (fFitList0Nt==0) {
        art::ServiceHandle<art::TFileService> tfs;
        fFitList0Nt = tfs->make<TNtuple>("fit0",
                                         "fitted track - det. frame",
                                         "run:subrun:evt:trk:x:y:z:slice");
      }
      nt = fFitList0Nt;
    }
    if (which==1) {
      if (fFitList1Nt==0) {
        art::ServiceHandle<art::TFileService> tfs;
        fFitList1Nt = tfs->make<TNtuple>("fit1",
                                         "fitted track - track frame",
                                         "run:subrun:evt:trk:x:y:z:slice");
      }
      nt = fFitList1Nt;
    }

    unsigned int i;
    float v[8];
    v[0] = fRun;
    v[1] = fSubrun;
    v[2] = fEvent;
    v[3] = fProng;
    v[7] = fSlice;
    for (i=0; i<z.size(); ++i) {
      v[4] = x[i];
      v[5] = y[i];
      v[6] = z[i];
      nt->Fill(v);
    }
  }
  
  //......................................................................
  //
  // Create an ntuple so we can examine the residuals. Variables are:
  // m   - measured value. x/y/alpha in the two views
  // sig - estimated 1-sigma error on measured value
  // fit - fitted value of parameter
  // w   - which kind of measurement is it. Code follows:
  //       0 = X measurement from x view
  //       1 = X measurement interpolated
  //       2 = Y measurement from the y view
  //       3 = Y measurement interpolated
  //       4 = angle "measurement" from the x view
  //       5 = angle "measurement" from the y view
  // zs  - for hits, the z value of the hit in the track frame,
  //       for angles, the s value along the track
  //
  void LogResiduals(bpfit::HitList3D&          hit,
                    bpfit::ScatteringSurfaces& scat,
                    bpfit::Lutz&               xfit,
                    bpfit::Lutz&               yfit,
                    double*                    alphax,
                    double*                    alphay,
                    float                      pdgid)
  {
    if (fResidNt==0) {
      art::ServiceHandle<art::TFileService> tfs;
      fResidNt = tfs->make<TNtuple>("resid", 
                                    "resid", 
                                    "run:subrun:evt:trk:m:sig:fit:w:pdgid:zs:slice");
    }
    unsigned int i;
    float v[11];
    v[0] = fRun;
    v[1] = fSubrun;
    v[2] = fEvent;
    v[3] = fProng;
    v[10]= fSlice;
    v[8] = pdgid;
    for (i=0; i<hit.size(); ++i) {
      v[4] = hit[i].fX;
      v[5] = hit[i].fSigX;
      v[6] = xfit.X(hit[i].fZ);
      v[9] = hit[i].fZ;
      if (hit[i].fView==geo::kX) v[7] = 0;
      if (hit[i].fView==geo::kY) v[7] = 1;
      fResidNt->Fill(v);

      v[4] = hit[i].fY;
      v[5] = hit[i].fSigY;
      v[6] = yfit.X(hit[i].fZ);
      v[9] = hit[i].fZ;
      if (hit[i].fView==geo::kY) v[7] = 2;
      if (hit[i].fView==geo::kX) v[7] = 3;
      fResidNt->Fill(v);
    }
    for (i=0; i<scat.N(); ++i) {
      v[4] = 0;
      v[5] = scat.SigAlpha(i);
      v[6] = alphax[i];
      v[7] = 4;
      v[9] = scat.S(i);
      fResidNt->Fill(v);

      v[4] = 0;
      v[5] = scat.SigAlpha(i);
      v[6] = alphay[i];
      v[7] = 5;
      v[9] = scat.S(i);
      fResidNt->Fill(v);
    }
 
  }

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

  void FitTracks(art::Ptr<rb::Vertex>& vertexp, 
                 art::Ptr<rb::Prong>&  prongp,
                 std::vector<rb::Track>& tracks,
                 std::vector<rb::FitSum>& fitsums)
  {
    unsigned int i, j;
    bpfit::TrackBasis sys;
    double p0[3], p1[3];

    bpfit::HitList3D hits0(prongp->NCell());
    bpfit::HitList3D hits1(prongp->NCell());
    std::vector<double> s(hits1.size());

    // This first try/catch block is to pick out hit list errors.  These
    // will cause all particle track fit assumptions to fail.
    try {
      // First check that this prong has hits in more than one plane in
      // each view. If not, we won't try to fit it.
      art::ServiceHandle<geo::Geometry> geom;    
      std::set<unsigned int> Xplanes;
      std::set<unsigned int> Yplanes;
      for(unsigned int p = 0; p < prongp->NCell(); ++p) {
        unsigned int plane = prongp->Cell(p).get()->Plane();
        if(geom->Plane(plane)->View() == geo::kX) Xplanes.insert(plane);
        if(geom->Plane(plane)->View() == geo::kY) Yplanes.insert(plane);
      }
      if(Xplanes.size() == 1 && Yplanes.size() == 1) {
        throw BPException(__FILE__, __LINE__, kNPLANES, 2);
      }
      if(Xplanes.size() == 1) {
        throw BPException(__FILE__, __LINE__, kNPLANES, 0);
      }
      if(Yplanes.size() == 1) {
        throw BPException(__FILE__, __LINE__, kNPLANES, 1);
      }

      //
      // Track fitting starts with a hit list
      //
      this->FillHitList(prongp, hits0);
      
      /*
      // Only average hits if either view has a small number of hits.
      if(prongp->NXCell()<10 || prongp->NYCell()<10){
        // Currently commented out the use of the averagexy function...
        // I have verified that for most fits nothing is affected but there are still
        // a few I need to dive deeper into.

        this->AverageXY(hits0);
      }
      */
    
      for(i=0; i<hits0.size(); ++i){
        hits1.SetHit(i,hits0[i].fView, hits0[i].fX, hits0[i].fY, hits0[i].fZ, hits0[i].fT,
                     hits0[i].fSigX, hits0[i].fSigY, hits0[i].fSigZ, hits0[i].fSigT,
                     hits0[i].fW);
      }
      
       if (fDebugHitList0) this->DebugHitList(hits0, 0);
      //
      // Then we want to swing around to a track-based coordinate system
      // so that we can properly locate the scattering planes
      //
      this->MakeBasis(vertexp, hits0, sys);

      for (i=0; i<hits0.size(); ++i) {
        sys.DetToTrk(hits0[i].fX,     hits0[i].fY,     hits0[i].fZ,
                     hits0[i].fSigX,  hits0[i].fSigY,  hits0[i].fSigZ,
                     &hits1[i].fX,    &hits1[i].fY,    &hits1[i].fZ,
                     &hits1[i].fSigX, &hits1[i].fSigY, &hits1[i].fSigZ);
      }
      if (fDebugHitList1) this->DebugHitList(hits1, 1);
      hits1.SortByZ();

      //
      // Get the extent of the track and z1 locations of measurement
      // planes
      //
      this->FindEndPoints(hits1, sys, p0, p1);

      for (i=0; i<hits1.size(); ++i) {
        s[i] = hits1[i].fZ;
      }
      
    } // end of try block
    catch (BPException& Ex) {
      if(fVerbose)      Ex.Print();
      if(fFillExceptNt) this->FillExceptNt(Ex);
      return;
    } // end of catch block

    //
    // Loop over the list of particle assumptions requested for this fit,
    // and make one track for each of those assumptions.
    //
    for(unsigned int t = 0; t < fPartList.size(); ++t) {
      // This second try/catch block is to pick out track fitting errors.
      // This will cause one particle fit assumption to fail but not all.
      try {

        // Make the track to be fit.
        //
        ///\ TODO:  For now we initialize it to the prong that it is fit from.
        //          We should be putting in the hits from the prongs with the
        //          weights used to calculate the track.
        rb::Track track(*prongp);
        rb::FitSum fitsum;

        art::ServiceHandle<geo::Geometry> geom;    
        bpfit::ScatteringSurfaces scat_surf(fPartList[t], fdStol, fdXtol, fdTtol);
        scat_surf.MakeSurfaces(*geom, p0, p1, s);

        if (fDebugScatSurf) this->DebugScatSurf(scat_surf, fPartList[t]);
        
        bpfit::Lutz lutzx(hits1.size(), scat_surf.N());
        bpfit::Lutz lutzy(hits1.size(), scat_surf.N());
        for (i=0; i<hits1.size(); ++i) {
          lutzx.SetMeasurement(i, hits1[i].fZ, hits1[i].fX, hits1[i].fSigX);
          lutzy.SetMeasurement(i, hits1[i].fZ, hits1[i].fY, hits1[i].fSigY);
        }
        for (j=0; j<scat_surf.N(); ++j) {
          lutzx.SetScatteringPlane(j, scat_surf.S(j), scat_surf.SigAlpha(j));
          lutzy.SetScatteringPlane(j, scat_surf.S(j), scat_surf.SigAlpha(j));
        }
        double alphax[scat_surf.N()];
        double alphay[scat_surf.N()];
        lutzx.Fit(0, 0, alphax, 0);
        lutzy.Fit(0, 0, alphay, 0);
        //
        // Check the fit results and remove far-away outliers
        //
        double chix, chiy, chim, fac;
        for (i=0; i<hits1.size(); ++i) {
          chix = lutzx.Chi2XIi(i);
          chiy = lutzy.Chi2XIi(i);
          chim = (chix>chiy) ? chix : chiy;
          fac  = sqrt(chim);
          if (chim>fChi2Cutoff) {
            hits1[i].fW = 1.0E-9;
            hits1[i].fSigX *= fac;
            hits1[i].fSigY *= fac;
            hits1[i].fSigZ *= fac;
          }
        }
        //
        // Need to recompute the projection from orthogonal projection as
        // outliers would have biased the estimates
        //
        hits1.CompleteOrtho();
        for (i=0; i<hits1.size(); ++i) {
          lutzx.SetMeasurement(i, hits1[i].fZ, hits1[i].fX, hits1[i].fSigX);
          lutzy.SetMeasurement(i, hits1[i].fZ, hits1[i].fY, hits1[i].fSigY);
        }
        lutzx.Fit(0, 0, alphax, 0);
        lutzy.Fit(0, 0, alphay, 0);
        
        //
        // Get the fitted track positions in the track frame. Take the
        // start and end points from the hits and insert the track
        // locations from the scattering surfaces.
        //
        std::vector<double> xfit(scat_surf.N()+2);
        std::vector<double> yfit(scat_surf.N()+2);
        std::vector<double> zfit(scat_surf.N()+2);
        zfit[0] = hits1[0].fZ;
        for (j=0; j<scat_surf.N(); ++j) {
          zfit[1+j] = scat_surf.S(j);
        }
        zfit[1+j] = hits1[hits1.size()-1].fZ;
        for (j=0; j<zfit.size(); ++j) {
          xfit[j] = lutzx.X(zfit[j]);
          yfit[j] = lutzy.X(zfit[j]);
        }
        
        if (fDebugFitList1) this->DebugFitList(1, xfit, yfit, zfit);
        if (fLogResiduals) {
          this->LogResiduals(hits1, scat_surf, lutzx, lutzy, alphax, alphay, fPartList[t]);
        }
        
        //
        // Translate back to the detector frame and build a path
        //
        std::vector<double> xfit_det;
        std::vector<double> yfit_det;
        std::vector<double> zfit_det;
        double xyz[3], dx, dy, dz;
        bpfit::Path path(*geom, zfit.size(), 0.0);
        for (j=0; j<zfit.size(); ++j) {
          sys.TrkToDet(xfit[j], yfit[j], zfit[j], 0,   0,   0,
                       &xyz[0], &xyz[1], &xyz[2], &dx, &dy, &dz);
          this->ClampToDetector(xyz);
          path.SetPoint(j, xyz);
          xfit_det.push_back(xyz[0]);
          yfit_det.push_back(xyz[1]);
          zfit_det.push_back(xyz[2]);
        }
        path.Scrub();
        if (fExtendUp!=0) path.AdjustStartPoint(fExtendUp);
        if (fExtendDn!=0) path.AdjustEndPoint(fExtendDn);

        // Debugging fit parameters in track frame.
        if (fDebugFitList0) this->DebugFitList(1, xfit, yfit, zfit);
        // Debugging fit parameters in detector frame.
        // if (fDebugFitList0) this->DebugFitList(0, xfit_det, yfit_det, zfit_det);
        
        //
        // Log the track parameters
        //
        if (fFillTrackNt) this->FillTrackNt(hits1, lutzx, lutzy, path, fPartList[t]);
        
        //
        // Build the track before returning
        //
        track.SetDir(path.Dir(0.0));
        track.SetStart(path.fXYZ[0]);   
        for (i=1; i<path.fXYZ.size(); ++i) {
          track.AppendTrajectoryPoint(path.fXYZ[i]);
        }
        
        // Calculate the degrees of freedom for each track.
        unsigned int nhitx = 0, nhity = 0;
        for (unsigned int i=0; i<hits1.size(); ++i) {
          if (hits1[i].fView==geo::kX && hits1[i].fW > 1.0E-9) ++nhitx;
          if (hits1[i].fView==geo::kY && hits1[i].fW > 1.0E-9) ++nhity;
        }

        // Check that the number of hits used in each view isn't too small.
        // The number of degrees of freedom for each track fit is N Hits
        // that went into making the track minus 4, so if used 0 hits in
        // either view or <= 4 total, we will throw this exception.
        if (nhitx < 1 || nhity < 1 || (nhitx+nhity) <= 4) {
          throw BPException(__FILE__, __LINE__, kBAD_NDOF, nhitx+nhity);
        }

        double momentum = 0.0;
        double mass     = 0.0;

        if(fPartList[t] == 13) {
          path.MuonParams(0, 0, &momentum, 0, 0);
          mass = 105.658E-3; 
        }
        if(fPartList[t] == 211) {
          path.PionParams(0, 0, &momentum, 0, 0);
          mass = 139.570E-3; 
        }
        if(fPartList[t] == 2212) {
          path.ProtonParams(0, 0, &momentum, 0, 0);
          mass = 938.272E-3;
        }

        // set the kinematic information and chi2 and ndof...
        fitsum.SetKin(path.Dir(0.0), momentum, mass);
        fitsum.SetPDG(fPartList[t]);
        fitsum.SetChi2(lutzx.Chi2XI(),lutzy.Chi2XI(),
                       lutzx.Chi2Beta(),lutzy.Chi2Beta(),
                       nhitx, nhity, scat_surf.N(), scat_surf.N());
        
        tracks.push_back(track);
        fitsums.push_back(fitsum);

      } // end of try block
      catch (BPException& Ex) {
        if(fVerbose)      Ex.Print();
        if(fFillExceptNt) this->FillExceptNt(Ex);
      } // end of catch block
    } // end loop over fPartList

    return;

  }

  //......................................................................
  
  void produce(art::Event& evt) 
  {
    //
    // Allocations for the track products and their associations
    //
    std::unique_ptr<std::vector<rb::Track> >
      trackcol(new std::vector<rb::Track>);
    std::unique_ptr<std::vector<rb::FitSum> >
      fitcol(new std::vector<rb::FitSum>);
    std::unique_ptr< art::Assns<rb::Track, rb::Prong> >
      trkPrngAssns(new art::Assns<rb::Track, rb::Prong>);
    std::unique_ptr< art::Assns<rb::Track, rb::Cluster> >
      trkSlcAssns(new art::Assns<rb::Track, rb::Cluster>);
    std::unique_ptr< art::Assns<rb::FitSum, rb::Track> >
      fitTrkAssns(new art::Assns<rb::FitSum, rb::Track>);

    fRun    = evt.run();
    fSubrun = evt.subRun();
    fEvent  = evt.id().event();
    fSlice  = 0;
    //
    // Extract the slices, vertices, and tracks we need as inputs
    //
    art::Handle< std::vector<rb::Cluster> > sliceh;
    art::PtrVector<rb::Cluster>             slicep;
    this->FindSlices(evt, sliceh, slicep);
    art::FindManyP<rb::Vertex> vertexh(sliceh, evt, fVertexLabel);
    art::FindManyP<rb::Prong>
      prongh(sliceh, evt, art::InputTag(fClusterLabel,fClusterInstance));
    if (!vertexh.isValid()) {
      return;
    }

    unsigned int islice;
    for (islice=0; islice<slicep.size(); ++islice) {
      fSlice = islice;

      //
      // Skip the noise slice - no physics there
      //
      if (slicep[islice]->IsNoise()) continue;

      //
      // To avoid problems with jobs running out of memory and/or taking too
      // long to process a single event, skip the slice if it has too many
      // hits since the reconstruction will be junk anyway.
      //
      if (slicep[islice]->NCell() > fMaxSliceHits) continue;

      //
      // Find vertex associated with this slice
      //
      std::vector< art::Ptr<rb::Vertex> > vertexp;
      vertexp = vertexh.at(islice);
      if (vertexp.size()!=1) continue;
      
      //
      // Find the tracks associated with this slice
      //
      std::vector< art::Ptr<rb::Prong> > prongp;
      this->FindProngs(evt, islice, prongh, vertexp, prongp);
      if (prongp.size()==0) continue;

      //
      // At this point we should have a list of tracks to fit. Go
      // ahead and process them
      //
      unsigned int nprong_good = 0; // How many tracks pass selection cuts
      for (fProng=0; fProng<prongp.size(); ++fProng) {
        //
        // Figure out if the track has enough information to be
        // fit. If not, skip to the next track
        //
        if (!this->GoodProng(prongp[fProng])) continue;
        ++nprong_good;

        //
        // Fit the tracks under the requested assumptions.
        //
        std::vector<rb::Track>   tracks;
        std::vector<rb::FitSum>  fitsums;

        this->FitTracks(vertexp[0], prongp[fProng], tracks, fitsums);
        
        // Add all returned tracks to the full list and make the appropriate Assns.
        assert(tracks.size() == fitsums.size());
        for(unsigned int t = 0; t < tracks.size(); ++t) {
          trackcol->push_back(tracks[t]);
          fitcol->push_back(fitsums[t]);
          util::CreateAssn(*this, evt, *trackcol, prongp[fProng], *trkPrngAssns);
          util::CreateAssn(*this, evt, *trackcol, slicep[islice] , *trkSlcAssns);
          util::CreateAssn(*this, evt, *fitTrkAssns, *fitcol, *trackcol);
        }    
                
      } // Loop on prongs
    } // Loop on slices

    evt.put(std::move(trackcol));
    evt.put(std::move(fitcol));
    evt.put(std::move(trkPrngAssns));
    evt.put(std::move(trkSlcAssns));
    evt.put(std::move(fitTrkAssns));

  } // produce() method
};

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

DEFINE_ART_MODULE(bpfit::BreakPoint)

////////////////////////////////////////////////////////////////////////