SplitTracks_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       SplitTracks
// Module Type: producer
// File:        SplitTracks_module.cc
//
// Generated at Tue Apr 23 11:01:33 2013 by Denver Whittington using artmod
// from art v1_02_06.
////////////////////////////////////////////////////////////////////////

// ART includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Optional/TFileService.h"

// C++ includes
#include <vector>
#include <string>
#include <cstdlib>
#include <iostream>

// NOvASoft includes
#include "CMap/service/CMapService.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Track.h"
#include "TrackFit/CosmicTrackAlg.h"
#include "TrackFit/WindowTrackingAlg.h"
#include "Utilities/AssociationUtil.h"

// ROOT includes
#include "TTree.h"
#include "TH1.h"
#include "TH2.h"

namespace align {
  class SplitTracks;
}

class align::SplitTracks : public art::EDAnalyzer {
public:
  explicit SplitTracks(fhicl::ParameterSet const &pset);
  virtual ~SplitTracks();

  void analyze(art::Event const &evt) override;
  //void produce(art::Event &evt);
  void beginJob() override;
  //void endJob() override;
  void ClearVectors();

  double DemingRegFit( const std::vector<double>& x,
                       const std::vector<double>& y,
                       const std::vector<double>& w,
                       double& sigXY, double& sigZ,
                       double& x1, double& y1,
                       double& x2, double& y2,
                       double& sigy1, double& sigy2 );

  TVector3 GeoEstimatePos( art::Ptr<rb::CellHit> cellHit, double W );
  void GetCellEndpoints( art::Ptr<rb::CellHit> cell, double *cellLo, double *cellHi );

private:

  std::string  fTrackInLabel;
  int          fMinBlocksInTrack;
  unsigned int fMinPlanesInSubCluster;

  trk::TrackAlg       *fTrkAlg;
  int                  fTrackAlgEnum;
  TTree *fGeoTree;
  TTree *fOutTree;

  // Geometry details
  std::vector<double> blockCenterX;
  std::vector<double> blockCenterY;
  std::vector<double> blockCenterZ;

  // Track details
  int trk_ID;
  double trk_StartX, trk_StartY, trk_StartZ;
  double trk_StopX, trk_StopY, trk_StopZ;
  std::vector<int> subtrk_Block;
  std::vector<double> subtrk_StartX, subtrk_StartY, subtrk_StartZ;
  std::vector<double> subtrk_ErrStartX, subtrk_ErrStartY, subtrk_ErrStartZ;
  std::vector<double> subtrk_StopX, subtrk_StopY, subtrk_StopZ;
  std::vector<double> subtrk_ErrStopX, subtrk_ErrStopY, subtrk_ErrStopZ;

  // Monitoring histograms
  TH1D *xzViewTrkRes;
  TH1D *yzViewTrkRes;
  TH2D *xzViewTrkResVsX;
  TH2D *xzViewTrkResVsZ;
  TH2D *yzViewTrkResVsY;
  TH2D *yzViewTrkResVsZ;
  TH1D *xzViewSubtrkRes;
  TH1D *yzViewSubtrkRes;
  TH2D *xzViewSubtrkResVsX;
  TH2D *xzViewSubtrkResVsZ;
  TH2D *yzViewSubtrkResVsY;
  TH2D *yzViewSubtrkResVsZ;
  TH1D *xzViewRefitRes;
  TH1D *yzViewRefitRes;
  TH2D *xzViewRefitResVsX;
  TH2D *xzViewRefitResVsZ;
  TH2D *yzViewRefitResVsY;
  TH2D *yzViewRefitResVsZ;

  TH1D *trackDCAtoCell;
  TH1D *subTrackDCAtoCell;
  TH1D *refitDCAtoCell;

  TH2D *subtrkXYstart;
  TH2D *subtrkXZstart;
  TH2D *subtrkYZstart;

};

// Constructor
align::SplitTracks::SplitTracks(fhicl::ParameterSet const &pset)
  : EDAnalyzer(pset)
  , fTrackInLabel(pset.get<std::string>("TrackInLabel"))
  , fMinBlocksInTrack(pset.get<int>("MinBlocksInTrack"))
  , fMinPlanesInSubCluster(pset.get<unsigned int>("MinPlanesInSubCluster"))
  , fTrkAlg(0)
  , fTrackAlgEnum(pset.get< int >("TrackingAlgEnumeration", 0))
{

  switch(fTrackAlgEnum){
  case 0:
    fTrkAlg = new trk::CosmicTrackAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
    break;
  case 1:
    fTrkAlg = new trk::WindowTrackingAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
    break;
  default:
    fTrkAlg = new trk::CosmicTrackAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
    break;
  }
}

// Destructor
align::SplitTracks::~SplitTracks()
{
  if (fTrkAlg) delete fTrkAlg;
}

// Begin job
void align::SplitTracks::beginJob() {
  art::ServiceHandle<art::TFileService> tfs;

  // Output TTree
  fGeoTree = tfs->make<TTree>("GeometryInfo","Geometry Info");
  fOutTree = tfs->make<TTree>("SplitTracks","Split Tracks");

  // Geometry details
  fGeoTree->Branch("blockCenterX", &blockCenterX);
  fGeoTree->Branch("blockCenterY", &blockCenterY);
  fGeoTree->Branch("blockCenterZ", &blockCenterZ);
  // Get Geo details...
  fGeoTree->Fill();

  // Track details
  fOutTree->Branch("trk_ID",    &trk_ID);
  fOutTree->Branch("trk_StartX", &trk_StartX);
  fOutTree->Branch("trk_StartY", &trk_StartY);
  fOutTree->Branch("trk_StartZ", &trk_StartZ);
  fOutTree->Branch("trk_StopX", &trk_StopX);
  fOutTree->Branch("trk_StopY", &trk_StopY);
  fOutTree->Branch("trk_StopZ", &trk_StopZ);

  fOutTree->Branch("subtrk_Block", &subtrk_Block);
  fOutTree->Branch("subtrk_StartX", &subtrk_StartX);
  fOutTree->Branch("subtrk_StartY", &subtrk_StartY);
  fOutTree->Branch("subtrk_StartZ", &subtrk_StartZ);
  fOutTree->Branch("subtrk_ErrStartX", &subtrk_ErrStartX);
  fOutTree->Branch("subtrk_ErrStartY", &subtrk_ErrStartY);
  fOutTree->Branch("subtrk_ErrStartZ", &subtrk_ErrStartZ);
  fOutTree->Branch("subtrk_StopX", &subtrk_StopX);
  fOutTree->Branch("subtrk_StopY", &subtrk_StopY);
  fOutTree->Branch("subtrk_StopZ", &subtrk_StopZ);
  fOutTree->Branch("subtrk_ErrStopX", &subtrk_ErrStopX);
  fOutTree->Branch("subtrk_ErrStopY", &subtrk_ErrStopY);
  fOutTree->Branch("subtrk_ErrStopZ", &subtrk_ErrStopZ);

  subtrkXYstart = tfs->make<TH2D>("subtrkXYstart","Downstream Subtrack Start (XY);X-view Cell;Y-view Cell",100,-775,775,100,-775,775);
  subtrkXZstart = tfs->make<TH2D>("subtrkXZstart","Downstream Subtrack Start (XZ);Plane;X-view Cell",600,0,6000,1550,-775,775);
  subtrkYZstart = tfs->make<TH2D>("subtrkYZstart","Downstream Subtrack Start (YZ);Plane;Y-view Cell",600,0,6000,1550,-775,775);

  xzViewTrkRes = tfs->make<TH1D>("xzViewTrkRes","Track Residual (X-View);Residual",100,-10,10);
  yzViewTrkRes = tfs->make<TH1D>("yzViewTrkRes","Track Residual (Y-View);Residual",100,-10,10);
  xzViewTrkResVsX = tfs->make<TH2D>("xzViewTrkResVsX","Track Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
  xzViewTrkResVsZ = tfs->make<TH2D>("xzViewTrkResVsZ","Track Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
  yzViewTrkResVsY = tfs->make<TH2D>("yzViewTrkResVsY","Track Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
  yzViewTrkResVsZ = tfs->make<TH2D>("yzViewTrkResVsZ","Track Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);

  xzViewSubtrkRes = tfs->make<TH1D>("xzViewSubtrkRes","SubTrack Residual (X-View);Residual",100,-10,10);
  yzViewSubtrkRes = tfs->make<TH1D>("yzViewSubtrkRes","SubTrack Residual (Y-View);Residual",100,-10,10);
  xzViewSubtrkResVsX = tfs->make<TH2D>("xzViewSubtrkResVsX","SubTrack Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
  xzViewSubtrkResVsZ = tfs->make<TH2D>("xzViewSubtrkResVsZ","SubTrack Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
  yzViewSubtrkResVsY = tfs->make<TH2D>("yzViewSubtrkResVsY","SubTrack Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
  yzViewSubtrkResVsZ = tfs->make<TH2D>("yzViewSubtrkResVsZ","SubTrack Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);

  xzViewRefitRes = tfs->make<TH1D>("xzViewRefitRes","Refitted SubTrack Residual (X-View);Residual",100,-10,10);
  yzViewRefitRes = tfs->make<TH1D>("yzViewRefitRes","Refitted SubTrack Residual (Y-View);Residual",100,-10,10);
  xzViewRefitResVsX = tfs->make<TH2D>("xzViewRefitResVsX","Refitted SubTrack Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
  xzViewRefitResVsZ = tfs->make<TH2D>("xzViewRefitResVsZ","Refitted SubTrack Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
  yzViewRefitResVsY = tfs->make<TH2D>("yzViewRefitResVsY","Refitted SubTrack Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
  yzViewRefitResVsZ = tfs->make<TH2D>("yzViewRefitResVsZ","Refitted SubTrack Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);

  trackDCAtoCell    = tfs->make<TH1D>("trackDCAtoCell","Track DCA to Cell",100,0,10);
  subTrackDCAtoCell = tfs->make<TH1D>("subTrackDCAtoCell","SubTrack DCA to Cell",100,0,10);
  refitDCAtoCell    = tfs->make<TH1D>("refitDCAtoCell","Refitted SubTrack DCA to Cell",100,0,10);

  return;
}

void align::SplitTracks::ClearVectors() {
  trk_ID = 0;
  trk_StartX = 0; trk_StartY = 0; trk_StartZ = 0;
  trk_StopX = 0; trk_StopY = 0; trk_StopZ = 0;
  subtrk_Block.clear();
  subtrk_StartX.clear();
  subtrk_StartY.clear();
  subtrk_StartZ.clear();
  subtrk_ErrStartX.clear();
  subtrk_ErrStartY.clear();
  subtrk_ErrStartZ.clear();
  subtrk_StopX.clear();
  subtrk_StopY.clear();
  subtrk_StopZ.clear();
  subtrk_ErrStopX.clear();
  subtrk_ErrStopY.clear();
  subtrk_ErrStopZ.clear();
  return;
}

// Produce
void align::SplitTracks::analyze(art::Event const &evt)
//void align::SplitTracks::produce(art::Event &evt)
{
  art::ServiceHandle<geo::Geometry> geoSvc;
  const int detId = geoSvc->DetId(); // What detector am I in?
  art::ServiceHandle<cmap::CMap> cmapSvc;

  // Get the tracks that were made by CosmicTrack
  // and split along alignment boundaries
  /* (blocks for now) */
  art::Handle< std::vector<rb::Track> > CosmicTracks;
  evt.getByLabel(fTrackInLabel,CosmicTracks);
  if ( CosmicTracks->empty() ) {
    std::cout << "No tracks found in " << fTrackInLabel << " for this event!" << std::endl;
    return;
  }

  // Loop through CosmicTracks
  std::unique_ptr< std::vector<rb::Track> > subTrackColl(new std::vector<rb::Track>);
  for ( size_t trkItr = 0; trkItr < CosmicTracks->size(); ++trkItr ) {
    art::Ptr<rb::Track> trackIn(CosmicTracks,trkItr);

    if ( !trackIn->Is3D() ) continue;
    this->ClearVectors();

    /* Track Requirements:
     *   Crosses a plane boundary
     *   Contains at least 3 hits in XZ and YZ in each plane
     */

    // Crosses a plane boundary?
    daqchannelmap::lchan logchan  = cmapSvc->Map()->encodeLChan(detId, trackIn->MaxPlane(), 1);
    daqchannelmap::lchan logchan1 = cmapSvc->Map()->encodeLChan(detId, trackIn->MinPlane(), 1); 
    int maxBlock = cmapSvc->Map()->computeBlock(logchan);
    int minBlock = cmapSvc->Map()->computeBlock(logchan1);
    if ( maxBlock < minBlock ) {
      std::cout << "MaxBlock is smaller than MinBlock for this track!" << std::endl;
      std::swap(maxBlock,minBlock);
      // This shouldn't happen, so throw a warning if it does...
    }
    if ( (maxBlock - minBlock + 1) < fMinBlocksInTrack ) continue;

    // Does the track go through the top and bottom of the detector bigbox?
    /*
    double start[3] = { trackIn->Start().X(), trackIn->Start().Y(), trackIn->Start().Z() };
    double dir[3]   = { trackIn->Dir().X(),   trackIn->Dir().Y(),   trackIn->Dir().Z()   };
    double boxMin[3];
    double boxMax[3];
    geoSvc->DetectorBigBox(&boxMin[0],&boxMax[0],&boxMin[1],&boxMax[1],&boxMin[2],&boxMax[2]);
    bool thruTop = geo::IntersectsBox( start, dir, boxMin[0], boxMax[0], boxMax[1], boxMax[1], boxMin[2], boxMax[2] );
    bool thruBot = geo::IntersectsBox( start, dir, boxMin[0], boxMax[0], boxMin[1], boxMin[1], boxMin[2], boxMax[2] );
    */
    //std::cout << "Min: " << trackIn->MinCell( geo::kY ) << "  Max: " << trackIn->MaxCell( geo::kY ) << std::endl;
    bool thruTop = ( trackIn->MaxCell( geo::kY ) == 32*12-1 );
    bool thruBot = ( trackIn->MinCell( geo::kY ) == 0     );
    if (!(thruTop&&thruBot) ) continue;

    // Record the parameters for this track
    trk_ID = trackIn->ID();
    trk_StartX = trackIn->Start().X();
    trk_StartY = trackIn->Start().Y();
    trk_StartZ = trackIn->Start().Z();
    trk_StopX = trackIn->Stop().X();
    trk_StopY = trackIn->Stop().Y();
    trk_StopZ = trackIn->Stop().Z();

    // Group CellHits into new clusters by block
    std::vector<rb::Cluster> subClusters( maxBlock-minBlock+1, rb::Cluster() );
    // Loop through hits in track and sort into subClusters by block
    for ( size_t cellItr = 0; cellItr < trackIn->NCell(); ++cellItr ) {
      const art::Ptr<rb::CellHit>& cell = trackIn->Cell(cellItr);
      daqchannelmap::lchan logchan2 = cmapSvc->Map()->encodeLChan(detId, cell->Plane(),1);
      int block = cmapSvc->Map()->computeBlock(logchan2);
      subClusters.at(block-minBlock).Add(cell); // WEIGHTS???
      /*
      double W = trackIn->W(cell.get());
      double cellPos[3];
      geoSvc->Plane(cell->Plane())->Cell(cell->Cell())->GetCenter(cellPos,W);
      */
      double cellLo[3], cellHi[3];
      this->GetCellEndpoints( cell, cellLo, cellHi );
      double dOnTrack, dInCell;
      TVector3 pcaOnTrack, pcaInCell;
      double dca = geo::ClosestApproach( trackIn->Start(), trackIn->Stop(), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
      trackDCAtoCell->Fill(dca);
      double cellPos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
      double trkX(0), trkY(0);
      trackIn->InterpolateXY(cellPos[2],trkX,trkY);
      if ( cell->View() == geo::kX ) {
        double res = trkX-cellPos[0];
        xzViewTrkRes->Fill(res);
        xzViewTrkResVsX->Fill(cell->Cell(),res);
        xzViewTrkResVsZ->Fill(cell->Plane(),res);
      }
      if ( cell->View() == geo::kY ) {
        double res = trkY-cellPos[1];
        yzViewTrkRes->Fill(res);
        yzViewTrkResVsY->Fill(cell->Cell(),res);
        yzViewTrkResVsZ->Fill(cell->Plane(),res);
      }
    }

    /*
    // Report!
    std::cout << "Hey! This track had " << trackIn->NCell() << " cell hits." << std::endl;
    std::cout << "     They were distributed across " << subClusters.size() << " blocks, with" << std::endl;
    for ( size_t i = 0; i < subClusters.size(); ++i ) {
      std::cout << "     " << subClusters[i].NCell() << " in block " << i << std::endl;
    }
    */

    // Loop through our subClusters and fit them to get end-points and errors
    for ( size_t clusItr = 0; clusItr < subClusters.size(); ++clusItr ) {

      // Make a track so we can refine its hit positions with the extrapolated W positions
      std::vector<rb::Track> subTracks = fTrkAlg->MakeTrack(&subClusters.at(clusItr));

      // only use subtracks if the vector is of size 1
      if(subTracks.size() < 1 || subTracks.size() > 1) continue;

      rb::Track subTrack = subTracks.front();

      // Does this cluster contain enough Planes?
      unsigned int nPlanes = subTrack.ExtentPlane();
      if ( nPlanes < fMinPlanesInSubCluster ) continue;
      if ( subTrack.NXCell() < 2 ) continue;
      if ( subTrack.NYCell() < 2 ) continue;
      subTrackColl->push_back(subTrack);

      // Refit the two views to get the uncertainties
      std::vector<double> xzViewX, xzViewZ, yzViewY, yzViewZ, xzViewW, yzViewW;
      double zMin = 10000000;
      double zMax = 0.;
      for ( size_t i = 0; i < subTrack.NCell(); ++i ) {
        art::Ptr<rb::CellHit> cellHit = subTrack.Cell(i);
        /*
        double W = subTrack.W(cellHit.get());
        double pos[3];
        geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->GetCenter(pos,W);
        */
        double cellLo[3], cellHi[3];
        this->GetCellEndpoints( cellHit, cellLo, cellHi );
        double dOnTrack, dInCell;
        TVector3 pcaOnTrack, pcaInCell;
        double dca = geo::ClosestApproach( subTrack.Start(), subTrack.Stop(), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
        subTrackDCAtoCell->Fill(dca);
        double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
        //double localXCheck[3]; geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->WorldToLocal(pos,localXCheck);
        //std::cout << localXCheck[0] << " " << localXCheck[1] << " " << localXCheck[2] << std::endl;
        double subtrkX(0), subtrkY(0);
        subTrack.InterpolateXY(pos[2],subtrkX,subtrkY);
        if ( cellHit->View() == geo::kX ) {
          double res = subtrkX-pos[0];
          xzViewSubtrkRes->Fill(res);
          xzViewSubtrkResVsX->Fill(cellHit->Cell(),res);
          xzViewSubtrkResVsZ->Fill(cellHit->Plane(),res);
          xzViewX.push_back(pos[0]);
          xzViewZ.push_back(pos[2]);
          xzViewW.push_back(1.);
        }
        if ( cellHit->View() == geo::kY ) {
          double res = subtrkY-pos[1];
          yzViewSubtrkRes->Fill(res);
          yzViewSubtrkResVsY->Fill(cellHit->Cell(),res);
          yzViewSubtrkResVsZ->Fill(cellHit->Plane(),res);
          yzViewY.push_back(pos[1]);
          yzViewZ.push_back(pos[2]);
          yzViewW.push_back(1.);
        }
        if ( pos[2] < zMin ) zMin = pos[2];
        if ( pos[2] > zMax ) zMax = pos[2];
      }

      // Fits
      double sigX = 3./sqrt(12);
      double sigY = 3./sqrt(12);
      double sigZ = 5./sqrt(12);
      double xStart, yStart, zStart;
      double xStartErr, yStartErr, zStartErr;
      double xStop, yStop, zStop;
      double xStopErr, yStopErr, zStopErr;
      zStart = zMin; zStartErr = 0.;//sigZ;
      zStop  = zMax; zStopErr  = 0.;//sigZ;
      this->DemingRegFit( xzViewZ, xzViewX, xzViewW, sigZ, sigX, zStart, xStart, zStop, xStop, xStartErr, xStopErr );
      this->DemingRegFit( yzViewZ, yzViewY, yzViewW, sigZ, sigY, zStart, yStart, zStop, yStop, yStartErr, yStopErr );

      // Residuals from refits
      for ( size_t i = 0; i < subTrack.NXCell(); ++i ) {
        art::Ptr<rb::CellHit> cellHit = subTrack.XCell(i);
        /*
        double W = subTrack.W(cellHit.get());
        double pos[3];
        geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->GetCenter(pos,W);
        */
        double cellLo[3], cellHi[3];
        this->GetCellEndpoints( cellHit, cellLo, cellHi );
        double dOnTrack, dInCell;
        TVector3 pcaOnTrack, pcaInCell;
        double dca = geo::ClosestApproach( TVector3(xStart,yStart,zStart), TVector3(xStop,yStop,zStop), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
        refitDCAtoCell->Fill(dca);
        double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
        double slopeX = (xStop-xStart)/(zStop-zStart);
        double extrapX = slopeX*(pos[2]-zStart)+xStart;
        double res = extrapX-pos[0];
        xzViewRefitRes->Fill(res);
        xzViewRefitResVsX->Fill(cellHit->Cell(),res);
        xzViewRefitResVsZ->Fill(cellHit->Plane(),res);
      }
      for ( size_t i = 0; i < subTrack.NYCell(); ++i ) {
        art::Ptr<rb::CellHit> cellHit = subTrack.YCell(i);
        /*
        double W = subTrack.W(cellHit.get());
        double pos[3];
        geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->GetCenter(pos,W);
        */
        double cellLo[3], cellHi[3];
        this->GetCellEndpoints( cellHit, cellLo, cellHi );
        double dOnTrack, dInCell;
        TVector3 pcaOnTrack, pcaInCell;
        double dca = geo::ClosestApproach( TVector3(xStart,yStart,zStart), TVector3(xStop,yStop,zStop), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
        refitDCAtoCell->Fill(dca);
        double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
        double slopeY = (yStop-yStart)/(zStop-zStart);
        double extrapY = slopeY*(pos[2]-zStart)+yStart;
        double res = extrapY-pos[1];
        yzViewRefitRes->Fill(res);
        yzViewRefitResVsY->Fill(cellHit->Cell(),res);
        yzViewRefitResVsZ->Fill(cellHit->Plane(),res);
      }
      subtrkXYstart->Fill(xStart,yStart);
      subtrkXZstart->Fill(zStart,xStart);
      subtrkYZstart->Fill(zStart,yStart);

      //if ( xStartErr > 1000 || xStopErr > 1000 ) continue;

      // Record the parameters for this subtrack
      daqchannelmap::lchan logchan3 = cmapSvc->Map()->encodeLChan(detId, subTrack.MinPlane(),1);
      subtrk_Block.push_back( cmapSvc->Map()->computeBlock(logchan3));
      subtrk_StartX.push_back( xStart );
      subtrk_StartY.push_back( yStart );
      subtrk_StartZ.push_back( zStart );
      subtrk_ErrStartX.push_back( xStartErr );
      subtrk_ErrStartY.push_back( yStartErr );
      subtrk_ErrStartZ.push_back( zStartErr );
      subtrk_StopX.push_back( xStop );
      subtrk_StopY.push_back( yStop );
      subtrk_StopZ.push_back( zStop );
      subtrk_ErrStopX.push_back( xStopErr );
      subtrk_ErrStopY.push_back( yStopErr );
      subtrk_ErrStopZ.push_back( zStopErr );

    }

    //std::cout << "This track had " << TracksOut.size() << " subtracks." << std::endl;

    // Save properties of this track and each of its subtracks to ROOT TTree
    fOutTree->Fill();

  } // End loop through CosmicTracks
  //evt.put(std::move(subTrackColl));

  return;
}

//......................................................................
///
/// \brief Find the best-fit line to a collection of points in 2-D
/// using Deming regression
///
/// This is similar to LinFitMinDPerp, but uses errors on hits
/// http://en.wikipedia.org/wiki/Deming_regression
///
/// @param x        - input vector of x coordinates
/// @param y        - input vector of y coordinates
/// @param w        - input vector of weights for the points
/// @param sigX     - uncertainty on cell hit position along x-axis
/// @param sigY     - uncertainty on cell hit position along y-axis
/// @param x1       - position of one end of segment (given)
/// @param x2       - position of other end of segment (given)
/// @param y1       - coordinate of line extrapolated to x1
/// @param y2       - coordinate of line extrapolated to x2
/// @param sigy1    - uncertainty on extrapolated coordinate y1
/// @param sigy2    - uncertainty on extrapolated coordinate y2
/// @returns The chi^2 value defined by chi^2 = sum_i[w_i d^2_i]
///
double align::SplitTracks::DemingRegFit(const std::vector<double>& x,
                                        const std::vector<double>& y,
                                        const std::vector<double>& w,
                                        double& sigX, double& sigY,
                                        double& x1, double& y1,
                                        double& x2, double& y2,
                                        double& sigy1, double& sigy2 )
{
  double delta = sigY / sigX;

  // Before going ahead, make sure we have sensible arrays
  assert(x.size()==y.size());
  assert(x.size()==w.size());
  assert(x.size()>=2);
  unsigned int i;

  // Accumulate the sums needed to compute the fit line
  double xAve = 0.0;
  double yAve = 0.0;
  int nHit = 0;
  for (i=0; i<w.size(); ++i) {
    if (w[i] == 0) continue;
    nHit++;
    xAve += x[i];
    yAve += y[i];
  }
  xAve /= double(nHit);
  yAve /= double(nHit);
  double Sxy = 0.0;
  double Syy = 0.0;
  double Sxx = 0.0;
  for (i=0; i<w.size(); ++i) {
    if (w[i] == 0) continue;
    Sxx += (x[i] - xAve)*(x[i] - xAve);
    Syy += (y[i] - yAve)*(y[i] - yAve);
    Sxy += (x[i] - xAve)*(y[i] - yAve);
  }
  Sxx /= double(nHit-1);
  Syy /= double(nHit-1);
  Sxy /= double(nHit-1);

  // The best fit through the points
  double M1;
  //double B1;
  double xydiff = Syy - delta*Sxx;
  double radic = sqrt( (xydiff*xydiff) + 4*delta*Sxy*Sxy );
  double numer = xydiff + radic;
  double denom = 2.*Sxy;
  M1 = ( numer ) / ( denom );
  //B1 = yAve - M1*xAve;

  y1 = M1*(x1-xAve) + yAve;
  y2 = M1*(x2-xAve) + yAve;

  // Compute errors
  double sig2y1(0), sig2y2(0);
  for (i=0; i<w.size(); ++i) {
    if (w[i] == 0) continue;
    double dnumdx = 0.;
    double dnumdy = 0.;
    dnumdx += delta*(2/double(nHit-1))*(x[i]-xAve);
    dnumdx += 1./(2*radic)*( 2*xydiff*( -delta*2/double(nHit-1)*(x[i]-xAve) )
                             + 8*delta*Sxy*1./double(nHit-1)*(y[i]-yAve) );
    dnumdy += (2/double(nHit-1))*(y[i]-yAve);
    dnumdx += 1./(2*radic)*( 2*xydiff*( 2/double(nHit-1)*(y[i]-yAve) )
                             + 8*Sxy*1./double(nHit-1)*(x[i]-xAve) );
    double dM1dx = denom != 0. ? ( dnumdx - M1*2/double(nHit-1)*(y[i]-yAve) ) / ( denom ) : 0.;
    double dM1dy = denom != 0. ? ( dnumdy - M1*2/double(nHit-1)*(x[i]-xAve) ) / ( denom ) : 0.;
    /*
    double dB1dx = -( M1/double(nHit) + xAve*dM1dx );
    double dB1dy = 1/double(nHit) - xAve*dM1dy;
    double dx1dx = 0;
    double dx2dx = 0;
    double dy1dx = dM1dx*x1 + M1*dx1dx + dB1dx;
    double dy1dy = dM1dy*x1 + dB1dy;
    double dy2dx = dM1dx*x2 + M1*dx2dx + dB1dx;
    double dy2dy = dM1dy*x2 + dB1dy;
    */
    double dy1dx = dM1dx*(x1-xAve) + M1;
    double dy1dy = dM1dy*(x1-xAve) + 1./double(nHit);
    double dy2dx = dM1dx*(x2-xAve) + M1;
    double dy2dy = dM1dy*(x2-xAve) + 1./double(nHit);
    sig2y1 += dy1dx*dy1dx*sigX*sigX + dy1dy*dy1dy*sigY*sigY;
    sig2y2 += dy2dx*dy2dx*sigX*sigX + dy2dy*dy2dy*sigY*sigY;
  }
  sigy1 = sqrt(sig2y1);
  sigy2 = sqrt(sig2y2);

  //if ( sigy1 > 1000 || sigy2 > 1000 ) {
  //  std::cout << "Anomalous solution! ( " << sigy1 << " " << sigy2 << " )" << std::endl;
  //  for ( i = 0; i < w.size(); ++i ) {
  //    if ( w[i] == 0 ) continue;
  //    std::cout << "(" << y[i] << "," << x[i] << ")  ";
  //  }
  //  std::cout << "Endpoints: (" << y1 << "," << x1 << ") (" << y2 << "," << x2 << ")" <<std::endl;
  //  std::cout << "Slope " << M1 << "  Intercept " << B1 << std::endl;
  //}

  // Compute the chi^2 function for return (Doesn't use errors!)
  double chi2 = 0.;

  return chi2;
}

TVector3 align::SplitTracks::GeoEstimatePos( art::Ptr<rb::CellHit> cellHit, double XY ) {
  art::ServiceHandle<geo::Geometry> geoSvc;
  size_t plane = cellHit->Plane();
  size_t cell = cellHit->Cell();
  const geo::CellGeo* theCell = geoSvc->Plane(plane)->Cell(cell);
  double localHitPos[3] = { 0., 0., XY };
  double worldHitPos[3];
  theCell->LocalToWorld(localHitPos,worldHitPos);
  return TVector3(worldHitPos);
}

void align::SplitTracks::GetCellEndpoints( art::Ptr<rb::CellHit> cellHit, double *cellLo, double *cellHi ) {
  art::ServiceHandle<geo::Geometry> geoSvc;
  size_t plane = cellHit->Plane();
  size_t cell = cellHit->Cell();
  const geo::CellGeo* theCell = geoSvc->Plane(plane)->Cell(cell);
  double localLo[3] = { 0., 0., -theCell->HalfL() };
  double localHi[3] = { 0., 0.,  theCell->HalfL() };
  theCell->LocalToWorld(localLo,cellLo);
  theCell->LocalToWorld(localHi,cellHi);
  return;
}

DEFINE_ART_MODULE(align::SplitTracks)