Track3DKalmanSPS_module.cc

Go to the documentation of this file.

//
// \file Track3DKalmanSPS.cxx
//
// \author echurch@fnal.gov
//
//  This algorithm is designed to reconstruct 3D tracks through
//  GENFIT's Kalman filter.

// C++ includes
#include <algorithm> // std::sort()
#include <cmath>
#include <iterator> // std::distance()
#include <sstream>
#include <string>
#include <vector>

// ROOT includes
#include "TDatabasePDG.h"
#include "TMath.h"
#include "TMatrixT.h"
#include "TPrincipal.h"
#include "TTree.h"
#include "TVector3.h"

// Framework includes
#include "fhiclcpp/ParameterSet.h"
#include "fhiclcpp/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"

// GENFIT includes
#include "larreco/Genfit/GFAbsTrackRep.h"
#include "larreco/Genfit/GFConstField.h"
#include "larreco/Genfit/GFFieldManager.h"
#include "larreco/Genfit/GFKalman.h"
#include "larreco/Genfit/GFTrack.h"
#include "larreco/Genfit/PointHit.h"
#include "larreco/Genfit/RKTrackRep.h"

// LArSoft includes
#include "larcore/Geometry/Geometry.h"
#include "larcore/Geometry/WireReadout.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Track.h"
//\todo Reconstruction Producers should never include SimulationBase objects
#include "lardata/Utilities/AssociationUtil.h"
#include "nusimdata/SimulationBase/MCTruth.h"

static bool sp_sort_3dz(const art::Ptr<recob::SpacePoint>& h1,
                        const art::Ptr<recob::SpacePoint>& h2)
{
  const double* xyz1 = h1->XYZ();
  const double* xyz2 = h2->XYZ();
  return xyz1[2] < xyz2[2];
}
static bool sp_sort_3dy(const art::Ptr<recob::SpacePoint>& h1,
                        const art::Ptr<recob::SpacePoint>& h2)
{
  const double* xyz1 = h1->XYZ();
  const double* xyz2 = h2->XYZ();
  return xyz1[1] < xyz2[1];
}
static bool sp_sort_3dx(const art::Ptr<recob::SpacePoint>& h1,
                        const art::Ptr<recob::SpacePoint>& h2)
{
  const double* xyz1 = h1->XYZ();
  const double* xyz2 = h2->XYZ();
  return xyz1[0] < xyz2[0];
}

static bool sp_sort_nsppts(const art::PtrVector<recob::SpacePoint>& h1,
                           const art::PtrVector<recob::SpacePoint>& h2)
{
  const unsigned int s1 = h1.size();
  const unsigned int s2 = h2.size();
  return s1 > s2;
}

namespace trkf {

  class Track3DKalmanSPS : public art::EDProducer {
  public:
    explicit Track3DKalmanSPS(fhicl::ParameterSet const& pset);

  private:
    void produce(art::Event& evt);
    void beginJob();
    void endJob();
    double energyLossBetheBloch(const double& mass, const double p);

    void rotationCov(TMatrixT<Double_t>& cov, const TVector3& u, const TVector3& v);
    std::vector<double> dQdxCalc(const art::FindManyP<recob::Hit>& h,
                                 const art::PtrVector<recob::SpacePoint>& s,
                                 const TVector3& p,
                                 const TVector3& d);

    std::string fClusterModuleLabel;  // label for input collection
    std::string fSpptModuleLabel;     // label for input collection
    std::string fGenieGenModuleLabel; // label for input MC single particle generator
    std::string fG4ModuleLabel;       // label for input MC single particle generator
    std::string fSortDim;             // direction in which to sort spacepoints

    //  TFile *fileGENFIT;
    TTree* tree;

    TMatrixT<Double_t>* stMCT;
    TMatrixT<Double_t>* covMCT;
    TMatrixT<Double_t>* stREC;
    TMatrixT<Double_t>* covREC;
    Float_t chi2;
    Float_t chi2ndf;
    int fcont;

    Float_t* fpRECL;
    Float_t* fpREC;
    Float_t* fpMCMom;
    Float_t* fpMCPos;
    Float_t* fState0;
    Float_t* fCov0;
    int nfail;
    int ndf;
    int nchi2rePass;
    int ispptvec;
    int nspptvec;
    unsigned int evtt;
    unsigned int nTrks;
    unsigned int fptsNo;
    Float_t* fshx;
    Float_t* fshy;
    Float_t* fshz;
    Float_t* feshx;
    Float_t* feshy;
    Float_t* feshz;
    Float_t* feshyz;
    Float_t* fupdate;
    Float_t* fchi2hit;
    Float_t* fth;
    Float_t* feth;
    Float_t* fedudw;
    Float_t* fedvdw;
    Float_t* feu;
    Float_t* fev;
    Float_t* fsep;
    Float_t* fdQdx;
    unsigned int fDimSize; // if necessary will get this from pset in constructor.
    Float_t* fPCmeans;
    Float_t* fPCevals;
    Float_t* fPCsigmas;
    Float_t* fPC1;
    Float_t* fPC2;
    Float_t* fPC3;

    std::vector<double> fPosErr;
    std::vector<double> fMomErr;
    std::vector<double> fMomStart;
    double fPerpLim;
    bool fDoFit;
    int fNumIt;
    uint16_t fMinNumSppts;
    double fErrScaleS;
    double fErrScaleM;
    int fDecimate;
    double fMaxUpdate;
    int fDecimateU;
    double fDistanceU;
    double fMaxUpdateU;
    double fMomLow;
    double fMomHigh;
    int fPdg;
    double fChi2Thresh;
    int fMaxPass;

    genf::GFAbsTrackRep* repMC;
    genf::GFAbsTrackRep* rep;

  }; // class Track3DKalmanSPS

}

namespace trkf {

  //-------------------------------------------------
  Track3DKalmanSPS::Track3DKalmanSPS(fhicl::ParameterSet const& pset)
    : EDProducer{pset}
    , fDoFit(true)
    , fNumIt(5)
    , fMinNumSppts(5)
    , fErrScaleS(1.0)
    , fErrScaleM(1.0)
    , fDecimate(1)
    , fMaxUpdate(0.10)
    , fDecimateU(1)
    , fDistanceU(1.)
    , fMaxUpdateU(0.10)
    , fMomLow(0.001)
    , fMomHigh(100.)
    , fPdg(-13)
    , fChi2Thresh(12.0E12)
    , fMaxPass(1)
  {
    fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
    fSpptModuleLabel = pset.get<std::string>("SpptModuleLabel");
    fGenieGenModuleLabel = pset.get<std::string>("GenieGenModuleLabel");
    fG4ModuleLabel = pset.get<std::string>("G4ModuleLabel");
    fPosErr = pset.get<std::vector<double>>("PosErr3");     // resolution. cm
    fMomErr = pset.get<std::vector<double>>("MomErr3");     // GeV
    fMomStart = pset.get<std::vector<double>>("MomStart3"); //
    fPerpLim = pset.get<double>("PerpLimit", 1.e6);         // PCA cut.
    fDoFit = pset.get<bool>("DoFit", true);                 // Der.
    fNumIt = pset.get<int>("NumIt", 5);                     // Number x2 passes per fit.
    fMinNumSppts =
      pset.get<int>("MinNumSppts", 5); // Min number of sppts in vector to bother fitting
    fErrScaleS = pset.get<double>("ErrScaleSim", 1.0);  // error scale.
    fErrScaleM = pset.get<double>("ErrScaleMeas", 1.0); // error scale.
    fDecimate = pset.get<int>("DecimateC", 40);         // Sparsify data.
    fMaxUpdate = pset.get<double>("MaxUpdateC", 0.1);   // 0-out.
    fDecimateU = pset.get<int>("DecimateU", 100);       // Sparsify data.
    fDistanceU =
      pset.get<double>("DistanceU", 10.0); // Require this separation on uncontained 2nd pass.
    fMaxUpdateU = pset.get<double>("MaxUpdateU", 0.02);       // 0-out.
    fMomLow = pset.get<double>("MomLow", 0.01);               // Fit Range.
    fMomHigh = pset.get<double>("MomHigh", 20.);              // Fit Range.
    fPdg = pset.get<int>("PdgCode", -13);                     // mu+ Hypothesis.
    fChi2Thresh = pset.get<double>("Chi2HitThresh", 12.0E12); //For Re-pass.
    fSortDim = pset.get<std::string>("SortDirection", "z");   // case sensitive
    fMaxPass = pset.get<int>("MaxPass", 2);                   // mu+ Hypothesis.
    bool fGenfPRINT;
    if (pset.get_if_present("GenfPRINT", fGenfPRINT)) {
      MF_LOG_WARNING("Track3DKalmanSPS_GenFit")
        << "Parameter 'GenfPRINT' has been deprecated.\n"
           "Please use the standard message facility to enable GenFit debug output.";
      // A way to enable debug output is all of the following:
      // - compile in debug mode (no optimization, no profiling)
      // - if that makes everything too noisy, add to have everything else quiet
      //   services.message.debugModules: [ "Track3DKalmanSPS" ]
      // - to print all the GenFit debug messages, set
      //   services.message.destinations.LogDebugFile.categories.Track3DKalmanSPS_GenFit.limit: -1
      //   (assuming there is a LogDebugFile destination already; for example
      //   see the settings in uboonecode/uboone/Utilities/services_microboone.fcl )
    }

    produces<std::vector<recob::Track>>();
    produces<art::Assns<recob::Track, recob::Cluster>>();
    produces<art::Assns<recob::Track, recob::SpacePoint>>();
    produces<art::Assns<recob::Track, recob::Hit>>();
  }

  //-------------------------------------------------
  // stolen, mostly, from GFMaterialEffects.
  double Track3DKalmanSPS::energyLossBetheBloch(const double& mass, const double p = 1.5)
  {
    const double charge(1.0);
    const double mEE(188.); // eV
    const double matZ(18.);
    const double matA(40.);
    const double matDensity(1.4);
    const double me(0.000511);

    double beta = p / std::sqrt(mass * mass + p * p);
    double gammaSquare = 1. / (1.0 - beta * beta);
    // 4pi.r_e^2.N.me = 0.307075, I think.
    double dedx = 0.307075 * matDensity * matZ / matA / (beta * beta) * charge * charge;
    double massRatio = me / mass;
    // me=0.000511 here is in GeV. So mEE comes in here in eV.
    double argument = gammaSquare * beta * beta * me * 1.E3 * 2. /
                      ((1.E-6 * mEE) * std::sqrt(1 + 2 * std::sqrt(gammaSquare) * massRatio +
                                                 massRatio * massRatio));

    if (mass == 0.0) return (0.0);
    if (argument <= exp(beta * beta)) { dedx = 0.; }
    else {
      dedx *= (log(argument) - beta * beta); // Bethe-Bloch [MeV/cm]
      dedx *= 1.E-3;                         // in GeV/cm, hence 1.e-3
      if (dedx < 0.) dedx = 0.;
    }
    return dedx;
  }

  void Track3DKalmanSPS::rotationCov(TMatrixT<Double_t>& cov, const TVector3& u, const TVector3& v)
  {
    TVector3 xhat(1.0, 0.0, 0.0);
    TVector3 yhat(0.0, 1.0, 0.0);
    TVector3 zhat(0.0, 0.0, 1.0);
    TVector3 w(u.Cross(v));
    TVector3 uprime(u);
    TVector3 vprime(w.Cross(xhat));  // vprime now lies in yz plane
    Double_t angle(v.Angle(vprime)); /* This is the angle through which v
                                       must rotate. */
    uprime.Rotate(angle, w);         // u now is rotated the same amount
    if (uprime * xhat < 0) {
      uprime.Rotate(TMath::Pi(), w);
      vprime.Rotate(TMath::Pi(), w);
      angle += TMath::Pi();
    }
    // Build the block-diagonal 5x5 matrix
    double c = TMath::Cos(angle), s = TMath::Sin(angle);
    TMatrixT<Double_t> rot(5, 5);
    rot[0][0] = 1.0;
    rot[1][1] = c;
    rot[1][2] = s;
    rot[2][1] = -s;
    rot[2][2] = c;
    rot[3][3] = c;
    rot[3][4] = s;
    rot[4][3] = -s;
    rot[4][4] = c;

    cov = rot * cov;
  }

  std::vector<double> Track3DKalmanSPS::dQdxCalc(const art::FindManyP<recob::Hit>& h,
                                                 const art::PtrVector<recob::SpacePoint>& s,
                                                 const TVector3& dir,
                                                 const TVector3& loc)
  {
    // For now just Collection plane.
    // We should loop over all views, more generally.
    geo::SigType_t sig(geo::kCollection);
    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
    static art::PtrVector<recob::SpacePoint>::const_iterator sstart(s.begin());
    art::PtrVector<recob::SpacePoint>::const_iterator sppt = s.begin();
    std::vector<double> v;

    double mindist(100.0); // cm
    auto spptminIt(sppt);
    while (sppt != s.end()) {
      if (((**sppt).XYZ() - loc).Mag() < mindist) {
        double dist = ((**sppt).XYZ() - loc).Mag();

        // Jump out if we're as close as 0.1 mm away.
        if (dist < mindist) {
          mindist = dist;
          spptminIt = sppt;
          if (mindist < 0.01) break;
        }
      }
      sppt++;
    }
    sstart = spptminIt; // for next time.
    unsigned int ind(std::distance(s.begin(), spptminIt));

    std::vector<art::Ptr<recob::Hit>> hitlist = h.at(ind);

    double wirePitch = 0.;
    double angleToVert = 0;
    double charge = 0.;

    // FIXME: Currently the last iteration sets the charge, wirePitch, and angleToVert variables.
    //        Is this expected?
    for (std::vector<art::Ptr<recob::Hit>>::const_iterator ihit = hitlist.begin();
         ihit != hitlist.end();
         ++ihit) {
      const recob::Hit& hit1 = **ihit;
      if (hit1.SignalType() != sig) continue;
      auto const& plane = wireReadoutGeom.Plane(hit1.WireID());
      charge = hit1.Integral();
      wirePitch = plane.WirePitch();
      angleToVert = plane.Wire(0).ThetaZ(false) - 0.5 * TMath::Pi();
    }

    double cosgamma =
      TMath::Abs(TMath::Sin(angleToVert) * dir.Y() + TMath::Cos(angleToVert) * dir.Z());

    v.push_back(charge / wirePitch / cosgamma);

    return v;
  }

  //-------------------------------------------------
  void Track3DKalmanSPS::beginJob()
  {

    art::ServiceHandle<art::TFileService const> tfs;

    stMCT = new TMatrixT<Double_t>(5, 1);
    covMCT = new TMatrixT<Double_t>(5, 5);
    stREC = new TMatrixT<Double_t>(5, 1);
    covREC = new TMatrixT<Double_t>(5, 5);

    fpMCMom = new Float_t[4];
    fpMCPos = new Float_t[4];
    fpREC = new Float_t[4];
    fpRECL = new Float_t[4];
    fState0 = new Float_t[5];
    fCov0 = new Float_t[25];
    fDimSize = 20000; // if necessary will get this from pset in constructor.

    fshx = new Float_t[fDimSize];
    fshy = new Float_t[fDimSize];
    fshz = new Float_t[fDimSize];
    feshx = new Float_t[fDimSize];
    feshy = new Float_t[fDimSize];
    feshz = new Float_t[fDimSize];
    feshyz = new Float_t[fDimSize];
    fupdate = new Float_t[fDimSize];
    fchi2hit = new Float_t[fDimSize];
    fth = new Float_t[fDimSize];
    feth = new Float_t[fDimSize];
    fedudw = new Float_t[fDimSize];
    fedvdw = new Float_t[fDimSize];
    feu = new Float_t[fDimSize];
    fev = new Float_t[fDimSize];
    fsep = new Float_t[fDimSize];
    fdQdx = new Float_t[fDimSize];

    fPC1 = new Float_t[3];
    fPC2 = new Float_t[3];
    fPC3 = new Float_t[3];
    fPCmeans = new Float_t[3];
    fPCsigmas = new Float_t[3];
    fPCevals = new Float_t[3];

    //   //TFile fileGENFIT("GENFITout.root","RECREATE");

    tree = tfs->make<TTree>("GENFITttree", "GENFITttree");
    //tree->Branch("stMCT",&stMCT,"stMCT[5]/F"); // "TMatrixT<Double_t>"

    //tree->Branch("stMCT","TMatrixD",&stMCT,64000,0);
    //tree->Branch("covMCT",covMCT,"covMCT[25]/F");
    tree->Branch("covMCT", "TMatrixD", &covMCT, 64000, 0);
    tree->Branch("stREC", fState0, "stREC[5]/F");
    //tree->Branch("stREC","TMatrixD",&stREC,64000,0);
    tree->Branch("covREC", fCov0, "covREC[25]/F");
    //tree->Branch("covREC","TMatrixD",&covREC,64000,0);

    tree->Branch("nchi2rePass", &nchi2rePass, "nchi2rePass/I");
    tree->Branch("ispptvec", &ispptvec, "ispptvec/I");
    tree->Branch("chi2", &chi2, "chi2/F");
    tree->Branch("nfail", &nfail, "nfail/I");
    tree->Branch("ndf", &ndf, "ndf/I");
    tree->Branch("evtNo", &evtt, "evtNo/I");
    tree->Branch("nspptvec", &nspptvec, "nspptvec/I");
    tree->Branch("chi2ndf", &chi2ndf, "chi2ndf/F");

    tree->Branch("trkNo", &nTrks, "trkNo/I");
    tree->Branch("ptsNo", &fptsNo, "ptsNo/I");
    tree->Branch("cont", &fcont, "cont/I"); //O? Yes, O. Not 0, not L, ...
    tree->Branch("shx", fshx, "shx[ptsNo]/F");
    tree->Branch("shy", fshy, "shy[ptsNo]/F");
    tree->Branch("shz", fshz, "shz[ptsNo]/F");
    tree->Branch("sep", fsep, "sep[ptsNo]/F");
    tree->Branch("dQdx", fdQdx, "dQdx[ptsNo]/F");
    tree->Branch("eshx", feshx, "eshx[ptsNo]/F");
    tree->Branch("eshy", feshy, "eshy[ptsNo]/F");
    tree->Branch("eshz", feshz, "eshz[ptsNo]/F");
    tree->Branch("eshyz", feshyz, "eshyz[ptsNo]/F");
    tree->Branch("update", fupdate, "update[ptsNo]/F");
    tree->Branch("chi2hit", fchi2hit, "chi2hit[ptsNo]/F");
    tree->Branch("th", fth, "th[ptsNo]/F");
    tree->Branch("eth", feth, "eth[ptsNo]/F");
    tree->Branch("edudw", fedudw, "edudw[ptsNo]/F");
    tree->Branch("edvdw", fedvdw, "edvdw[ptsNo]/F");
    tree->Branch("eu", feu, "eu[ptsNo]/F");
    tree->Branch("ev", fev, "ev[ptsNo]/F");

    tree->Branch("pcMeans", fPCmeans, "pcMeans[3]/F");
    tree->Branch("pcSigmas", fPCsigmas, "pcSigmas[3]/F");
    tree->Branch("pcEvals", fPCevals, "pcEvals[3]/F");
    tree->Branch("pcEvec1", fPC1, "pcEvec1[3]/F");
    tree->Branch("pcEvec2", fPC2, "pcEvec2[3]/F");
    tree->Branch("pcEvec3", fPC3, "pcEvec3[3]/F");

    tree->Branch("pMCMom", fpMCMom, "pMCMom[4]/F");
    tree->Branch("pMCPos", fpMCPos, "pMCPos[4]/F");
    tree->Branch("pRECKalF", fpREC, "pRECKalF[4]/F");
    tree->Branch("pRECKalL", fpRECL, "pRECKalL[4]/F");
  }

  //-------------------------------------------------
  void Track3DKalmanSPS::endJob()
  {
    if (!rep) delete rep;
    if (!repMC) delete repMC;

    delete[] fpREC;
    delete[] fpRECL;
    delete[] fState0;
    delete[] fCov0;

    delete[] fshx;
    delete[] fshy;
    delete[] fshz;
    delete[] feshx;
    delete[] feshy;
    delete[] feshyz;
    delete[] feshz;
    delete[] fupdate;
    delete[] fchi2hit;
    delete[] fth;
    delete[] feth;
    delete[] fedudw;
    delete[] fedvdw;
    delete[] feu;
    delete[] fev;
    delete[] fsep;
    delete[] fdQdx;

    delete[] fPCmeans;
    delete[] fPCsigmas;
    delete[] fPCevals;
    delete[] fPC1;
    delete[] fPC2;
    delete[] fPC3;
  }

  //------------------------------------------------------------------------------------//
  void Track3DKalmanSPS::produce(art::Event& evt)
  {
    rep = 0;
    repMC = 0;
    // get services
    art::ServiceHandle<geo::Geometry const> geom;

    // Make a std::unique_ptr<> for the thing you want to put into the event
    // because that handles the memory management for you
    std::unique_ptr<std::vector<recob::Track>> tcol(new std::vector<recob::Track>);
    std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> tspassn(
      new art::Assns<recob::Track, recob::SpacePoint>);
    std::unique_ptr<art::Assns<recob::Track, recob::Hit>> thassn(
      new art::Assns<recob::Track, recob::Hit>);

    unsigned int tcnt = 0;

    // get input Hit object(s).
    art::Handle<std::vector<recob::Cluster>> clusterListHandle;
    evt.getByLabel(fClusterModuleLabel, clusterListHandle);

    art::Handle<std::vector<art::PtrVector<recob::SpacePoint>>> spptListHandle;
    evt.getByLabel(fSpptModuleLabel, spptListHandle);

    art::PtrVector<simb::MCTruth> mclist;

    if (!evt.isRealData()) {

      art::Handle<std::vector<simb::MCTruth>> mctruthListHandle;
      evt.getByLabel(fGenieGenModuleLabel, mctruthListHandle);

      for (unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii) {
        art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle, ii);
        mclist.push_back(mctparticle);
      }
    }

    std::vector<art::PtrVector<recob::SpacePoint>> spptIn(spptListHandle->begin(),
                                                          spptListHandle->end());
    // Get the spptvectors that are largest to be first, and smallest last.
    std::sort(spptIn.begin(), spptIn.end(), sp_sort_nsppts);

    std::vector<art::PtrVector<recob::SpacePoint>>::const_iterator sppt = spptIn.begin();
    auto spptB = sppt;

    TVector3 MCOrigin;
    TVector3 MCMomentum;
    // TVector3 posErr(.05,.05,.05); // resolution. 0.5mm
    // TVector3 momErr(.1,.1,0.2);   // GeV
    TVector3 posErr(fPosErr[0], fPosErr[1], fPosErr[2]);    // resolution. 0.5mm
    TVector3 momErr(fMomErr[0], fMomErr[1], fMomErr[2]);    // GeV
    TVector3 momErrFit(fMomErr[0], fMomErr[1], fMomErr[2]); // GeV

    // This is strictly for MC
    if (!evt.isRealData()) {
      // Below breaks are stupid, I realize. But rather than keep all the MC
      // particles I just take the first primary, e.g., muon and only keep its
      // info in the branches of the Ttree. I could generalize later, ...
      for (unsigned int ii = 0; ii < mclist.size(); ++ii) {
        art::Ptr<simb::MCTruth> mc(mclist[ii]);
        for (int jj = 0; jj < mc->NParticles(); ++jj) {
          simb::MCParticle part(mc->GetParticle(jj));
          MCOrigin.SetXYZ(part.Vx(), part.Vy(), part.Vz()); // V for Vertex
          MCMomentum.SetXYZ(part.Px(), part.Py(), part.Pz());
          MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
            << "FROM MC TRUTH, the particle's pdg code is: " << part.PdgCode()
            << " with energy = " << part.E() << ", with energy = " << part.E()
            << "\n  vtx: " << genf::ROOTobjectToString(MCOrigin)
            << "\n  momentum: " << genf::ROOTobjectToString(MCMomentum)
            << "\n    (both in Global (not volTPC) coords)";

          repMC = new genf::RKTrackRep(MCOrigin, MCMomentum, posErr, momErr, part.PdgCode());
          break;
        }
        break;
      }
      //for saving of MC truth
      stMCT->ResizeTo(repMC->getState());
      *stMCT = repMC->getState();
      covMCT->ResizeTo(repMC->getCov());
      *covMCT = repMC->getCov();
      MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
        << " repMC, covMC are ... \n"
        << genf::ROOTobjectToString(repMC->getState()) << genf::ROOTobjectToString(repMC->getCov());

    } // !isRealData
    nTrks = 0;
    TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(fPdg);
    Double_t mass = part->Mass();

    auto const& tpc = geom->TPC({0, 0});
    size_t cntp(0);
    while (sppt != spptIn.end()) {

      const art::PtrVector<recob::SpacePoint>& spacepoints = *sppt;

      double fMaxUpdateHere(fMaxUpdateU);
      int fDecimateHere(fDecimateU);
      double fErrScaleSHere(fErrScaleS);
      double fErrScaleMHere(fErrScaleM);
      int rePass0(1);
      unsigned int nTailPoints = 0; // 100;
      if (spacepoints.size() < 5) {
        sppt++;
        rePass0 = 3;
        continue;
      } // for now...

      MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
        << "\n\t found " << spacepoints.size()
        << " 3D spacepoint(s) for this element of std::vector<art:PtrVector> spacepoints. \n";

      // Let's find track's principle components.
      // We will sort along that direction, rather than z.
      // Further, we will skip outliers away from main axis.

      TPrincipal* principal = new TPrincipal(3, "ND");

      // I need to shuffle these around, so use copy constructor
      // to make non-const version spacepointss.
      art::PtrVector<recob::SpacePoint> spacepointss(spacepoints);

      // What I need is a nearest neighbor sorting.
      if (fSortDim.compare("y") && fSortDim.compare("x"))
        std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dz);
      if (!fSortDim.compare("y")) std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dy);
      if (!fSortDim.compare("x")) std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dx);

      for (unsigned int point = 0; point < spacepointss.size(); ++point) {
        if (point < (spacepointss.size() - nTailPoints)) {
          principal->AddRow(spacepointss[point]->XYZ());
        }
      }
      principal->MakePrincipals();

      const TVectorD* evals = principal->GetEigenValues();
      const TMatrixD* evecs = principal->GetEigenVectors();
      const TVectorD* means = principal->GetMeanValues();
      const TVectorD* sigmas = principal->GetSigmas();

      Double_t tmp[3], tmp2[3];
      principal->X2P((Double_t*)(means->GetMatrixArray()), tmp);
      principal->X2P((Double_t*)(sigmas->GetMatrixArray()), tmp2);
      for (unsigned int ii = 0; ii < 3; ++ii) {
        fPCmeans[ii] = (Float_t)(tmp[ii]);
        fPCsigmas[ii] = (Float_t)(tmp2[ii]);
        fPCevals[ii] = (Float_t)(evals->GetMatrixArray())[ii];
        // This method requires apparently pulling all 9
        // elements. Maybe 3 works.
        // Certainly, w can't be a scalar, I discovered.
        double w[9];
        evecs->ExtractRow(ii, 0, w);
        fPC1[ii] = w[0];
        fPC2[ii] = w[1];
        fPC3[ii] = w[2];
      }
      Double_t tmp3[3], tmp4[3], tmp5[3];
      principal->X2P((Double_t*)fPC1, tmp3);
      principal->X2P((Double_t*)fPC2, tmp4);
      principal->X2P((Double_t*)fPC3, tmp5);

      // Use a mip approximation assuming straight lines
      // and a small angle wrt beam.
      fMomStart[0] = spacepointss[spacepointss.size() - 1]->XYZ()[0] - spacepointss[0]->XYZ()[0];
      fMomStart[1] = spacepointss[spacepointss.size() - 1]->XYZ()[1] - spacepointss[0]->XYZ()[1];
      fMomStart[2] = spacepointss[spacepointss.size() - 1]->XYZ()[2] - spacepointss[0]->XYZ()[2];
      // This presumes a 0.8 GeV/c particle
      double dEdx = energyLossBetheBloch(mass, 1.0);
      // mom is really KE.
      TVector3 mom(dEdx * fMomStart[0], dEdx * fMomStart[1], dEdx * fMomStart[2]);
      double pmag2 = pow(mom.Mag() + mass, 2. - mass * mass);
      mom.SetMag(std::sqrt(pmag2));
      // Over-estimate by just enough for contained particles (5%).
      mom.SetMag(1.0 * mom.Mag());
      // If 1st/last point is close to edge of TPC, this track is
      // uncontained.Give higher momentum starting value in
      // that case.
      bool uncontained(false);
      double close(5.);     // cm.
      double epsMag(0.001); // cm.
      double epsX(250.0);   // cm.
      double epsZ(0.001);   // cm.

      if (spacepointss[spacepointss.size() - 1]->XYZ()[0] > (2. * tpc.HalfWidth() - close) ||
          spacepointss[spacepointss.size() - 1]->XYZ()[0] < close ||
          spacepointss[0]->XYZ()[0] > (2. * tpc.HalfWidth() - close) ||
          spacepointss[0]->XYZ()[0] < close ||
          spacepointss[spacepointss.size() - 1]->XYZ()[1] > (1. * tpc.HalfHeight() - close) ||
          (spacepointss[spacepointss.size() - 1]->XYZ()[1] < -1. * tpc.HalfHeight() + close) ||
          spacepointss[0]->XYZ()[1] > (1. * tpc.HalfHeight() - close) ||
          spacepointss[0]->XYZ()[1] < (-1. * tpc.HalfHeight() + close) ||
          spacepointss[spacepointss.size() - 1]->XYZ()[2] > (tpc.Length() - close) ||
          spacepointss[spacepointss.size() - 1]->XYZ()[2] < close ||
          spacepointss[0]->XYZ()[2] > (tpc.Length() - close) || spacepointss[0]->XYZ()[2] < close)
        uncontained = true;
      fErrScaleSHere = fErrScaleS;
      fErrScaleMHere = fErrScaleM;

      if (uncontained) {
        // Big enough to not run out of gas right at end of
        // track and give large angular deviations which
        // will kill the fit.
        mom.SetMag(2.0 * mom.Mag());
        MF_LOG_DEBUG("Track3DKalmanSPS_GenFit") << "Uncontained track ... ";
        fDecimateHere = fDecimateU;
        fMaxUpdateHere = fMaxUpdateU;
      }
      else {
        MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
          << "Contained track ... Run " << evt.run() << " Event " << evt.id().event();
        // Don't decimate contained tracks as drastically,
        // and omit only very large corrections ...
        // which hurt only high momentum tracks.
        fDecimateHere = fDecimate;
        fMaxUpdateHere = fMaxUpdate;
      }
      fcont = (int)(!uncontained);

      // This seems like best place to jump back to for a re-pass.
      unsigned short rePass = rePass0; // 1 by default;
      unsigned short maxPass(fMaxPass);
      unsigned short tcnt1(0);
      while (rePass <= maxPass) {

        TVector3 momM(mom);
        TVector3 momErrFit(momM[0] / 3.0, momM[1] / 3.0,
                           momM[2] / 3.0); // GeV

        genf::GFFieldManager::getInstance()->init(new genf::GFConstField(0.0, 0.0, 0.0));
        genf::GFDetPlane planeG((TVector3)(spacepointss[0]->XYZ()), momM);

        // Initialize with 1st spacepoint location and ...
        rep = new genf::RKTrackRep((TVector3)(spacepointss[0]->XYZ()),
                                   momM,
                                   posErr,
                                   momErrFit,
                                   fPdg); // mu+ hypothesis

        genf::GFTrack fitTrack(rep); //initialized with smeared rep
        fitTrack.setPDG(fPdg);
        // Gonna sort in z cuz I want to essentially transform here to volTPC coords.
        // volTPC coords, cuz that's what the Geant3/Geane stepper wants, as that's its understanding
        // from the Geant4 geometry, which it'll use. EC, 7-Jan-2011.
        int ihit = 0;
        fptsNo = 0;
        std::vector<unsigned int> spptSurvivedIndex;
        std::vector<unsigned int> spptSkippedIndex;
        unsigned int ppoint(0);
        for (unsigned int point = 0; point < spacepointss.size(); ++point) {
          double sep;
          // Calculate the distance in 2nd and 3rd PCs and
          // reject spt if it's too far out. Remember, the
          // sigmas are std::sqrt(eigenvals).
          double tmp[3];
          principal->X2P((Double_t*)(spacepointss[point]->XYZ()), tmp);
          sep = std::sqrt(tmp[1] * tmp[1] / fPCevals[1] + tmp[2] * tmp[2] / fPCevals[2]);
          if ((std::abs(sep) > fPerpLim) && (point < (spacepointss.size() - nTailPoints)) &&
              rePass <= 1) {
            spptSkippedIndex.push_back(point);
            continue;
          }
          // If point is too close in Mag or Z or too far in X from last kept point drop it.
          // I think this is largely redundant with PCA cut.
          TVector3 one(spacepointss[point]->XYZ());
          TVector3 two(spacepointss[ppoint]->XYZ());
          if (rePass == 2 && uncontained) {
            epsMag = fDistanceU; // cm
            fNumIt = 4;
            fErrScaleMHere = 0.1;
            // Above allows us to pretend as though measurements
            // are perfect, which we can ostensibly do now with
            // clean set of sppts. This creates larger gains, bigger
            // updates: bigger sensitivity to multiple scattering.
          }
          else if (rePass == 2 && !uncontained) {}
          if (point > 0 &&
              ((one - two).Mag() < epsMag ||               // too close
               ((one - two).Mag() > 8.0 && rePass == 1) || // too far
               std::abs(spacepointss[point]->XYZ()[2] - spacepointss[ppoint]->XYZ()[2]) < epsZ ||
               std::abs(spacepointss[point]->XYZ()[0] - spacepointss[ppoint]->XYZ()[0]) > epsX)) {
            spptSkippedIndex.push_back(point);
            continue;
          }

          if (
            point % fDecimateHere &&
            rePass <=
              1) // Jump out of loop except on every fDecimate^th pt. fDecimate==1 never sees continue.
          {
            /* Replace continue with a counter that will be used to index into vector of GFKalman fits. */
            continue;
          }

          ppoint = point;
          TVector3 spt3 = (TVector3)(spacepointss[point]->XYZ());
          std::vector<double> err3;
          err3.push_back(spacepointss[point]->ErrXYZ()[0]);
          err3.push_back(spacepointss[point]->ErrXYZ()[2]);
          err3.push_back(spacepointss[point]->ErrXYZ()[4]);
          err3.push_back(spacepointss[point]->ErrXYZ()[5]); // lower triangle diags.
          if (fptsNo < fDimSize) {
            fshx[fptsNo] = spt3[0];
            fshy[fptsNo] = spt3[1];
            fshz[fptsNo] = spt3[2];
            feshx[fptsNo] = err3[0];
            feshy[fptsNo] = err3[1];
            feshz[fptsNo] = err3[3];
            feshyz[fptsNo] = err3[2];
            fsep[fptsNo] = sep;

            if (fptsNo > 1) {
              TVector3 pointer(fshx[fptsNo] - fshx[fptsNo - 1],
                               fshy[fptsNo] - fshy[fptsNo - 1],
                               fshz[fptsNo] - fshz[fptsNo - 1]);
              TVector3 pointerPrev(fshx[fptsNo - 1] - fshx[fptsNo - 2],
                                   fshy[fptsNo - 1] - fshy[fptsNo - 2],
                                   fshz[fptsNo - 1] - fshz[fptsNo - 2]);
              fth[fptsNo] = (pointer.Unit()).Angle(pointerPrev.Unit());
            }
            feth[fptsNo] = 0.0;
            fedudw[fptsNo] = 0.0;
            fedvdw[fptsNo] = 0.0;
            feu[fptsNo] = 0.0;
            fev[fptsNo] = 0.0;
            fupdate[fptsNo] = 0.0;
          }

          MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
            << "ihit xyz..." << spt3[0] << "," << spt3[1] << "," << spt3[2];

          fitTrack.addHit(new genf::PointHit(spt3, err3),
                          1, //dummy detector id
                          ihit++);
          spptSurvivedIndex.push_back(point);
          fptsNo++;
        } // end loop over spacepoints.

        if (fptsNo <=
            fMinNumSppts) // Cuz 1st 2 in each direction don't count. Should have, say, 3 more.
        {
          MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
            << "Bailing cuz only " << fptsNo << " spacepoints.";
          rePass++;
          continue;
        }
        MF_LOG_DEBUG("Track3DKalmanSPS_GenFit") << "Fitting on " << fptsNo << " spacepoints.";

        genf::GFKalman k;
        k.setBlowUpFactor(5);   // 500 out of box. EC, 6-Jan-2011.
        k.setMomHigh(fMomHigh); // Don't fit above this many GeV.
        k.setMomLow(fMomLow);   // Don't fit below this many GeV.

        k.setInitialDirection(+1); // Instead of 1 out of box. EC, 6-Jan-2011.
        k.setNumIterations(fNumIt);
        k.setMaxUpdate(fMaxUpdateHere); // 0 out abs(update) bigger than this.
        k.setErrorScaleSTh(fErrScaleSHere);
        k.setErrorScaleMTh(fErrScaleMHere);

        bool skipFill = false;
        std::vector<TMatrixT<double>> hitMeasCov;
        std::vector<TMatrixT<double>> hitUpdate;
        std::vector<TMatrixT<double>> hitCov;
        std::vector<TMatrixT<double>> hitCov7x7;
        std::vector<TMatrixT<double>> hitState;
        std::vector<double> hitChi2;
        std::vector<TVector3> hitPlaneXYZ;
        std::vector<TVector3> hitPlaneUxUyUz;
        std::vector<TVector3> hitPlaneU;
        std::vector<TVector3> hitPlaneV;

        try {
          if (fDoFit) k.processTrack(&fitTrack);
        }
        catch (cet::exception& e) {
          MF_LOG_ERROR("Track3DKalmanSPS")
            << "just caught a cet::exception: " << e.what()
            << "\nExceptions won't be further handled; skip filling big chunks of the TTree.";
          skipFill = true;
        }

        if (rep->getStatusFlag() == 0) // 0 is successful completion
        {
          if (mf::isDebugEnabled()) {

            std::ostringstream dbgmsg;
            dbgmsg << "Original plane:";
            planeG.Print(dbgmsg);

            dbgmsg << "Current (fit) reference Plane:";
            rep->getReferencePlane().Print(dbgmsg);

            dbgmsg << "Last reference Plane:";
            rep->getLastPlane().Print(dbgmsg);

            if (planeG != rep->getReferencePlane())
              dbgmsg << "  => original hit plane (not surprisingly) not current reference Plane!";

            MF_LOG_DEBUG("Track3DKalmanSPS_GenFit") << dbgmsg.str();
          }
          if (!skipFill) {
            hitMeasCov = fitTrack.getHitMeasuredCov();
            hitUpdate = fitTrack.getHitUpdate();
            hitCov = fitTrack.getHitCov();
            hitCov7x7 = fitTrack.getHitCov7x7();
            hitState = fitTrack.getHitState();
            hitChi2 = fitTrack.getHitChi2();
            hitPlaneXYZ = fitTrack.getHitPlaneXYZ();
            hitPlaneUxUyUz = fitTrack.getHitPlaneUxUyUz();
            hitPlaneU = fitTrack.getHitPlaneU();
            hitPlaneV = fitTrack.getHitPlaneV();
            unsigned int totHits = hitState.size();

            // Pick up info from last fwd Kalman pass.
            unsigned int jhit = 0;
            for (unsigned int ihit = totHits - 2 * totHits / (2 * fNumIt);
                 ihit < (totHits - totHits / (2 * fNumIt));
                 ihit++) // was ihit<ihit<(totHits-fptsNo)<7
            {
              feth[jhit] = (Float_t)(hitMeasCov.at(ihit)[0][0]); // eth
              fedudw[jhit] = (Float_t)(hitMeasCov.at(ihit)[1][1]);
              fedvdw[jhit] = (Float_t)(hitMeasCov.at(ihit)[2][2]);
              feu[jhit] = (Float_t)(hitMeasCov.at(ihit)[3][3]);
              fev[jhit] = (Float_t)(hitMeasCov.at(ihit)[4][4]);
              fupdate[jhit] = (Float_t)(hitUpdate.at(ihit)[0][0]);
              fchi2hit[jhit] = (Float_t)(hitChi2.at(ihit));
              jhit++;
            }

            stREC->ResizeTo(rep->getState());
            *stREC = rep->getState();
            covREC->ResizeTo(rep->getCov());
            *covREC = rep->getCov();
            double dum[5];
            double dum2[5];
            for (unsigned int ii = 0; ii < 5; ii++) {
              stREC->ExtractRow(ii, 0, dum);
              fState0[ii] = dum[0];
              covREC->ExtractRow(ii, 0, dum2);
              for (unsigned int jj = 0; jj < 5; jj++) {
                fCov0[ii * 5 + jj] = dum2[jj];
              }
            }
            MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
              << " First State and Cov:" << genf::ROOTobjectToString(*stREC)
              << genf::ROOTobjectToString(*covREC);
            chi2 = (Float_t)(rep->getChiSqu());
            ndf = rep->getNDF();
            nfail = fitTrack.getFailedHits();
            nchi2rePass = (int)rePass;
            ispptvec = 1 + std::distance(spptB, sppt);
            chi2ndf = (Float_t)(chi2 / ndf);

            nTrks++;
            MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
              << "Track3DKalmanSPS about to do tree->Fill(). Chi2/ndf is " << chi2 / ndf << ".";
            fpMCMom[3] = MCMomentum.Mag();
            for (int ii = 0; ii < 3; ++ii) {
              fpMCMom[ii] = MCMomentum[ii];
              fpMCPos[ii] = MCOrigin[ii];
              fpREC[ii] = hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 * fNumIt))[ii];
              fpRECL[ii] = hitPlaneUxUyUz.at(totHits - totHits / (2 * fNumIt) - 1)[ii];
            }

            evtt = (unsigned int)evt.id().event();
            nspptvec = (unsigned int)spptListHandle->size();

            cntp++;
            std::vector<std::vector<double>> dQdx;
            // Calculate LastFwdPass quantities.
            std::vector<TMatrixT<double>> hitCovLFP;
            std::vector<TVector3> hitPlaneXYZLFP;
            std::vector<TVector3> hitPlaneUxUyUzLFP;
            std::vector<TVector3> hitPlanePxPyPzLFP;
            std::vector<TVector3> hitPlaneULFP;
            std::vector<TVector3> hitPlaneVLFP;
            std::vector<double> pLFP;
            std::vector<TMatrixT<double>> c7x7LFP;

            art::FindManyP<recob::Hit> hitAssns(spacepointss, evt, fSpptModuleLabel);
            for (unsigned int ii = 0; ii < totHits / (2 * fNumIt); ii++) {
              pLFP.push_back(1. / hitState.at(totHits - 2 * totHits / (2 * fNumIt) + ii)[0][0]);
              // hitCov -> hitCov7x7 !! EC, 11-May-2012.
              c7x7LFP.push_back(hitCov7x7.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
              hitCovLFP.push_back(hitCov.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
              hitPlaneXYZLFP.push_back(hitPlaneXYZ.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
              hitPlaneUxUyUzLFP.push_back(
                hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
              hitPlanePxPyPzLFP.push_back(hitPlaneUxUyUzLFP.back() * pLFP.back());
              hitPlaneULFP.push_back(hitPlaneU.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
              hitPlaneVLFP.push_back(hitPlaneV.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
              // Transform cov appropriate for track rotated
              // about w, forcing
              // v to be in y-z plane and u pointing in
              // +-ive x direction, per TrackAna convention.

              rotationCov(hitCovLFP.back(), hitPlaneULFP.back(), hitPlaneVLFP.back());
              dQdx.push_back(
                dQdxCalc(hitAssns, spacepointss, hitPlaneUxUyUzLFP.back(), hitPlaneXYZLFP.back()));
              fdQdx[ii] = dQdx.back().back();
            }
            fpREC[3] = rep->getMom(rep->getReferencePlane()).Mag();
            fpRECL[3] = pLFP[1];

            tree->Fill();

            // Put newest track on stack for this set of sppts,
            // remove previous one.
            recob::tracking::SMatrixSym55 covVtx, covEnd;
            for (unsigned int i = 0; i < 5; i++) {
              for (unsigned int j = i; j < 5; j++) {
                covVtx(i, j) = hitCovLFP.front()(i, j);
                covEnd(i, j) = hitCovLFP.back()(i, j);
              }
            }
            recob::Track the3DTrack(
              recob::TrackTrajectory(recob::tracking::convertCollToPoint(hitPlaneXYZLFP),
                                     recob::tracking::convertCollToVector(hitPlanePxPyPzLFP),
                                     recob::Track::Flags_t(hitPlaneXYZLFP.size()),
                                     true),
              0,
              -1.,
              0,
              covVtx,
              covEnd,
              tcnt++);
            if (rePass == 1) tcnt1++; // won't get here if Trackfit failed.
            if (rePass != 1 && tcnt1) tcol->pop_back();
            tcol->push_back(the3DTrack);
            util::CreateAssn(evt, *tcol, spacepointss, *tspassn);
            art::PtrVector<recob::Hit> hits; // = hitAssns;
            for (unsigned int ii = 0; ii < spacepointss.size(); ++ii) {
              hits.insert(hits.end(), hitAssns.at(ii).begin(), hitAssns.at(ii).end());
            }
            util::CreateAssn(evt, *tcol, hits, *thassn, tcol->size() - 1);
          } // end !skipFill
        }   // getStatusFlag

        if (!rep) delete rep;
        rePass++;
        // need to first excise bad spacepoints.
        // Grab up large Chi2hits first.
        art::PtrVector<recob::SpacePoint> spacepointssExcise;
        for (unsigned int ind = 0; ind < spptSurvivedIndex.size(); ++ind) {
          // Stricter to chuck sppts from uncontained then contained trks.
          if ((uncontained && fchi2hit[ind] > fChi2Thresh) ||
              (!uncontained && fchi2hit[ind] > 1.e9) || fchi2hit[ind] < 0.0
              // =0 eliminates ruled-out large updates. Not obviously
              // helpful.
              // add a restriction on dQdx here ...
          ) {
            art::PtrVector<recob::SpacePoint>::iterator spptIt =
              spacepointss.begin() + spptSurvivedIndex[ind];
            spacepointssExcise.push_back(*spptIt);
          }
        }
        // Now grab up those sppts which we skipped and don't want
        // to reconsider cuz they're too close to each other or
        // cuz they're too far in x, e.g.
        for (unsigned int ind = 0; ind < spptSkippedIndex.size(); ++ind) {
          art::PtrVector<recob::SpacePoint>::iterator spptIt =
            spacepointss.begin() + spptSkippedIndex[ind];
          spacepointssExcise.push_back(*spptIt);
        }
        // Get rid of redundantly Excised sppts before proceeding.
        std::stable_sort(spacepointss.begin(), spacepointss.end());
        std::stable_sort(spacepointssExcise.begin(), spacepointssExcise.end());
        std::set_union(spacepointssExcise.begin(),
                       spacepointssExcise.end(),
                       spacepointssExcise.begin(),
                       spacepointssExcise.end(),
                       spacepointssExcise.begin());
        // Now excise. New spacepointss will be smaller for second pass.
        art::PtrVector<recob::SpacePoint>::iterator diffSpptIt =
          std::set_difference(spacepointss.begin(),
                              spacepointss.end(),
                              spacepointssExcise.begin(),
                              spacepointssExcise.end(),
                              spacepointss.begin());
        spacepointss.erase(diffSpptIt, spacepointss.end());

        // calculate new seed momentum, and errors as merited
        if (rePass == 2 /* && uncontained */) {
          if (fpREC[3] < fMomHigh && fpREC[3] > fMomLow) {
            double kick(0.9); //Try to get away with a smaller start
            // for contained tracks. While for uncontained tracks
            // let's start up at a higher momentum and come down.
            // Though, 2 (1) GeV/c tracks are too low (high), so
            // instead let's actually lower starting value on
            // this second pass. -- EC 7-Mar-2013
            if (uncontained) kick = 0.5;
            for (int ii = 0; ii < 3; ++ii) {
              mom[ii] = momM[ii] * kick;
            }
          }
          else if (uncontained) {
            double unstick(1.0);
            if (fpREC[3] >= fMomHigh) unstick = 0.3;
            for (int ii = 0; ii < 3; ++ii) {
              mom[ii] = momM[ii] * unstick;
            }
          }
          else
            for (int ii = 0; ii < 3; ++ii) {
              mom[ii] = 1.1 * momM[ii];
            }
        }

      } // end while rePass<=maxPass

      sppt++;

    } // end while on elements of std::vector of art::PtrVectors of SpPts.

    if (!repMC) delete repMC;

    evt.put(std::move(tcol));
    // and now the spacepoints
    evt.put(std::move(tspassn));
    // and the hits. Note that these are all the hits from all the spacepoints considered,
    // even though they're not all contributing to the tracks.
    evt.put(std::move(thassn));
  }

  DEFINE_ART_MODULE(Track3DKalmanSPS)

} // end namespace