CosmicTrackAlg.cxx

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
/// \file CosmicTrackAlg.cxx
/// \brief   Track finder for cosmic rays
/// \author  brebel@fnal.gov, gsdavies@iastate.edu
/// \date
///////////////////////////////////////////////////////////////////////
#include "TrackFit/CosmicTrackAlg.h"

// NOvASoft includes
#include "RecoBase/Track.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "Utilities/func/MathUtil.h"

// ART includes
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

namespace trk {

  // Constructor
  CosmicTrackAlg::CosmicTrackAlg(fhicl::ParameterSet const &pset)
    : TrackAlg(pset)
  {
    fDHitGood = pset.get<double>("DHitGood");
  }

  // Destructor
  CosmicTrackAlg::~CosmicTrackAlg() {}

  // MakeTrack
  std::vector<rb::Track> CosmicTrackAlg::MakeTrack( const rb::Cluster* slice ) {

    std::vector<rb::Track> tracks;

    // We'll need some information about the
    // cell locations. That comes from the Geometry object.
    //
    art::ServiceHandle<geo::Geometry> geom;

    // flags for failure modes
    bool flagxz = true; 
    bool flagyz = true;

    // loop over all the hits for this slice and get the
    // x or y positions for the views

    std::vector<double> xzViewX;
    std::vector<double> xzViewZ;
    std::vector<double> yzViewY;
    std::vector<double> yzViewZ;
    std::vector<double> xzViewW; // W = weight
    std::vector<double> yzViewW; // W = weight
    double xyz[3]; // xyz[0]=x, xyz[1]=y, xyz[2]=z
    const int hitMax = slice->NCell();
    for(int hitIdx = 0; hitIdx < hitMax; ++hitIdx){
      art::Ptr<rb::CellHit> chit = slice->Cell(hitIdx);

      double sig = 1.;
      geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
      if(chit->View() == geo::kX){
        xzViewX.push_back(xyz[0]);
        xzViewZ.push_back(xyz[2]);
        xzViewW.push_back(sig);

        LOG_DEBUG("CosmicTrack") 
          << "x view " << xzViewZ[xzViewZ.size()-1] 
          << " " << xzViewX[xzViewX.size()-1] 
          << " " << chit->TNS();
      }
      if(chit->View() == geo::kY){
        yzViewY.push_back(xyz[1]);
        yzViewZ.push_back(xyz[2]);
        yzViewW.push_back(sig);
          
        LOG_DEBUG("CosmicTrack") 
          << "y view " << yzViewZ[yzViewZ.size()-1] 
          << " " << yzViewY[yzViewY.size()-1] 
          << " " << chit->TNS();
      }
    }// end loop over hits in this track
      
    if(xzViewZ.size() < 2 && yzViewZ.size() < 2){
      LOG_DEBUG("CosmicTrack") 
        << " Not enough hits in either view to fit for track";
      return tracks;
    }

    // Want to use the geo::LinFitMinDperp method to figure which hits are
    // closest to the line.  Need vectors of x/z, y/z and uncertainties for
    // the positions and it outputs the start and end points for the line in
    // each 2D view.  Function returns the chi^2
        
    double xzChi2        = 0.;
    double xzStart[2]    = {0.}; 
    double xzEnd[2]      = {0.}; 
    unsigned int xzNDrop = this->FitView(xzViewZ, xzViewX, xzViewW, xzChi2, xzStart, xzEnd);
    if(xzEnd[0]==xzStart[0] && xzEnd[1]==xzStart[1] && xzViewW.size()>1) flagxz = false ; 
    if(xzViewZ.size() == 1){
      xzStart[0] = xzEnd[0] = xzViewZ[0];
      xzStart[1] = xzEnd[1] = xzViewX[0];
    }

    double yzChi2        = 0.;
    double yzStart[2]    = {0.}; 
    double yzEnd[2]      = {0.}; 
    unsigned int yzNDrop = this->FitView(yzViewZ, yzViewY, yzViewW, yzChi2, yzStart, yzEnd);
    if(yzEnd[0]==yzStart[0] && yzEnd[1]==yzStart[1] && yzViewW.size()>1) flagyz = false ; 
    if(yzViewZ.size() == 1){
      yzStart[0] = yzEnd[0] = yzViewZ[0];
      yzStart[1] = yzEnd[1] = yzViewY[0];
    }

    LOG_DEBUG("CosmicTrack") 
      << "# xz hits = " << xzViewZ.size()-xzNDrop 
      << " chi^2 = " << xzChi2 << std::endl
      << " # yz hits = " << yzViewZ.size()-yzNDrop 
      << " chi^2 = " << yzChi2;


    // Tracks with only plane in a view are going to give infinite gradients
    // when we try to extend them later. Tilt them over as much as possible.
    // Try to make line extrapolate into detector rather than towards the outside.
    // Especially important for case of fits with no change in z in either view.

    if(fabs(xzEnd[0]-xzStart[0]) < 1){ // check xz view for 'verticalness'

      double ZaveY = yzEnd[0]-yzStart[0]/2.; // Center position of Z in yz view

      if( fabs(xzStart[1]) < fabs(xzEnd[1]) ) { // make sure to lean 'into' detector
        if( xzStart[0] < ZaveY ) { // make sure to point it the right direction
          xzStart[0] += geom->Plane(0)->Cell(0)->HalfD();
          xzEnd[0]   -= geom->Plane(0)->Cell(0)->HalfD(); }
        if( xzStart[0] > ZaveY ) { 
          xzStart[0] -= geom->Plane(0)->Cell(0)->HalfD();
          xzEnd[0]   += geom->Plane(0)->Cell(0)->HalfD(); }
      }
      else { // other end is closer to center
        if( xzEnd[0] > ZaveY ) { 
          xzStart[0] += geom->Plane(0)->Cell(0)->HalfD();
          xzEnd[0]   -= geom->Plane(0)->Cell(0)->HalfD(); }
        if( xzEnd[0] < ZaveY ) { 
          xzStart[0] -= geom->Plane(0)->Cell(0)->HalfD();
          xzEnd[0]   += geom->Plane(0)->Cell(0)->HalfD(); }
      }
    }

    if(fabs(yzEnd[0]-yzStart[0]) < 1){ // check yz view for verticalness

      double ZaveX = xzEnd[0]-xzStart[0]/2.; // Center position of Z in xz view

      if( fabs(yzStart[1]) < fabs(yzEnd[1]) ) { // make sure to lean 'into' detector
        if( yzStart[0] < ZaveX ) { // make sure to point it the right direction
          yzStart[0] += geom->Plane(0)->Cell(0)->HalfD();
          yzEnd[0]   -= geom->Plane(0)->Cell(0)->HalfD(); }
        if( yzStart[0] > ZaveX ) { 
          yzStart[0] -= geom->Plane(0)->Cell(0)->HalfD();
          yzEnd[0]   += geom->Plane(0)->Cell(0)->HalfD(); }
      }
      else { // other end is closer to center
        if( yzEnd[0] > ZaveX ) { 
          yzStart[0] += geom->Plane(0)->Cell(0)->HalfD();
          yzEnd[0]   -= geom->Plane(0)->Cell(0)->HalfD(); }
        if( yzEnd[0] < ZaveX ) { 
          yzStart[0] -= geom->Plane(0)->Cell(0)->HalfD();
          yzEnd[0]   += geom->Plane(0)->Cell(0)->HalfD(); }
      }
    }    

    // which z values define the end points?
    double startz = TMath::Min(yzStart[0], xzStart[0]);
    double endz   = TMath::Max(yzEnd[0],   xzEnd[0]);

    // are there at least 2 hits in this view? if not, reset the start
    // and end z positions to be based on the view with hits
    if(xzViewZ.size() < 2){
      //        startz = yzStart[0];
      //        endz   = yzEnd[0];
      flagxz=false;
    }
    if(yzViewZ.size() < 2){
      //        startz = xzStart[0];
      //        endz   = xzEnd[0];
      flagyz=false;
    }

    if(flagxz&&flagyz) {

      // Adjust the start and end points of each view to reflect how far the
      // z extent got moved
      const double xzSlope = (xzEnd[1]-xzStart[1])/(xzEnd[0]-xzStart[0]);
      xzStart[1] += (startz-xzStart[0])*xzSlope;
      xzEnd[1] += (endz-xzEnd[0])*xzSlope;

      const double yzSlope = (yzEnd[1]-yzStart[1])/(yzEnd[0]-yzStart[0]);
      yzStart[1] += (startz-yzStart[0])*yzSlope;
      yzEnd[1] += (endz-yzEnd[0])*yzSlope;

      // Now it's safe to update the z endpoints themselves
      xzStart[0] = yzStart[0] = startz;
      xzEnd[0] = yzEnd[0] = endz;

      LOG_DEBUG("CosmicTrack") 
        << startz << " " << endz << std::endl
        << xzStart[0] << " " << xzStart[1] << " " 
        << xzEnd[0] << " " << xzEnd[1] << std::endl
        << yzStart[0] << " " << yzStart[1] << " " 
        << yzEnd[0] << " " << yzEnd[1] << std::endl;

      // swap the start and end if timing indicates the track goes upstream
      if(!IsTrackDownstreamFromTiming(slice)){
        std::swap(xzStart[0], xzEnd[0]);
        std::swap(xzStart[1], xzEnd[1]);
        std::swap(yzStart[0], yzEnd[0]);
        std::swap(yzStart[1], yzEnd[1]);
        std::swap(startz, endz);
      }


      // now go through each view and only keep hits that lie near the line in
      // that view (ie FitView hasn't set their weight to zero).
      art::PtrVector<rb::CellHit> chits;

      //find min and max planes
      int minPlane = 9999;
      int maxPlane = -9999;

      for(unsigned int hitIdx = 0; hitIdx < slice->NXCell(); ++hitIdx){
        if(xzViewW[hitIdx] != 0){
          art::Ptr<rb::CellHit> ch = slice->XCell(hitIdx);
          chits.push_back(ch);
          if (ch->Plane() < minPlane) minPlane = ch->Plane();
          if (ch->Plane() > maxPlane) maxPlane = ch->Plane();
        }
      }
      for(unsigned int hitIdx = 0; hitIdx < slice->NYCell(); ++hitIdx){
        if(yzViewW[hitIdx] != 0){
          art::Ptr<rb::CellHit> ch = slice->YCell(hitIdx);
          chits.push_back(ch);
          if (ch->Plane() < minPlane) minPlane = ch->Plane();
          if (ch->Plane() > maxPlane) maxPlane = ch->Plane();
        }
      }

      if(chits.empty()){
        throw art::Exception(art::errors::ProductNotFound)
          << "No hits left on the track after clean up stage ";
      }

      //sort cells by plane, low to high before doing trajectory calculation
      rb::SortByPlane(chits);
      //if track is starting from the upstream end, reverse the sort
      if (startz > endz) std::reverse(chits.begin(),chits.end());
      
      // We can write out a proper 3D track; assume downward going
      TVector3 start(xzStart[1], yzStart[1], startz);
      TVector3 stop(xzEnd[1], yzEnd[1], endz); 

      if(start[1] < stop[1]) {
        TVector3 temp = start;
        start = stop;
        stop = temp;      
      }

      //sort cells by plane, low to high before doing trajectory calculation
      rb::SortByPlane(chits);
      //if track is starting from the upstream end, reverse the sort
      if (start[2] > stop[2]) std::reverse(chits.begin(),chits.end());

      //restrict track to exist inside detector; if this isn't possible, the track is no good, so don't write
      if(geo::ClampRayToDetectorVolume(&start, &stop, &(*geom))) { 

        rb::Track track(chits, start, stop-start);

        //find valid z range incase ends of tracks got clipped
        double minZ = std::min(start[2],stop[2]);
        double maxZ = std::max(start[2],stop[2]);

        // add a trajectory point for each plane - set
        // the point based on the z position of each plane 
        // and the direction of the track
        // the first and last planes are the start and stop
        // points, presumably

        int prevPlane = chits.front()->Plane();
        double xyz[3] = {0.};
        for(size_t c = 1; c < chits.size()-1; ++c){
          int plane = chits.at(c)->Plane();

          //check in case there are more then one hit in min/max plane not to add point multiple times
          if ((plane == minPlane) || (plane == maxPlane)) continue;

          if(plane != prevPlane) {
            geom->Plane(plane)->Cell(chits.at(c)->Cell())->GetCenter(xyz);
            //make sure plane is within range of track if ends were clipped to the detector volume
            if ((xyz[2] > minZ) && (xyz[2] < maxZ)){
              TVector3 dir = track.Dir().Unit();
              TVector3 point(start.X() + (xyz[2] - start.Z())*dir.X()/dir.Z(), 
                             start.Y() + (xyz[2] - start.Z())*dir.Y()/dir.Z(),
                             xyz[2]);
              track.AppendTrajectoryPoint(point);
            }
          }

          prevPlane = plane;
        }
        
        track.AppendTrajectoryPoint(stop);

        tracks.push_back(track);
        return tracks;
      }
      // if track isn't good in this case, just don't write it out      

    }

    // if at least one of the views was bad, no 3-D track.  Write out any 2-D tracks that are worthwhile
    if(flagxz&&!flagyz) { //xz view track only

      art::PtrVector<rb::CellHit> chits;

      for(unsigned int hitIdx = 0; hitIdx < slice->NXCell(); ++hitIdx){
        if(xzViewW[hitIdx] != 0) chits.push_back(slice->XCell(hitIdx));
      }

      rb::Track track(chits, geo::kX, xzStart[1], xzStart[0], xzEnd[1]-xzStart[1], xzEnd[0]-xzStart[0]);
      track.AppendTrajectoryPoint(xzEnd[1],xzEnd[0]);
      tracks.push_back(track);
      return tracks;
    }

    if(flagyz&&!flagxz) { //yz view track only
      art::PtrVector<rb::CellHit> chits;

      for(unsigned int hitIdx = 0; hitIdx < slice->NYCell(); ++hitIdx){
        if(yzViewW[hitIdx] != 0) chits.push_back(slice->YCell(hitIdx));
      }

      rb::Track track(chits, geo::kY, yzStart[1], yzStart[0], yzEnd[1]-yzStart[1], yzEnd[0]-yzStart[0]);
      track.AppendTrajectoryPoint(yzEnd[1],yzEnd[0]);
      tracks.push_back(track);
      return tracks;
    }

    return tracks;
  } // End of MakeTrack

  //......................................................................
  unsigned int CosmicTrackAlg::FitView(std::vector<double> &x,
                                       std::vector<double> &y,
                                       std::vector<double> &w_orig,
                                       double &chiSqr,
                                       double *start,
                                       double *end)
  {
    // We need to store the original weights in case we need to restore them
    std::vector<double> w = w_orig;

    const unsigned int nhit = w.size();

    chiSqr = -9999.;

    // First pass, drop hits that are clear outliers (more than 2x tolerance
    // from anything else). If it turns out they really are good hits (eg
    // surrounded by bad channels) they will be added back in by the last
    // step. But this helps prevent distant noise hits from throwing off the
    // original fit.
    unsigned int ndrop = 0;
    for(unsigned int i = 0; i < nhit; ++i){
      double closest = 1e10;
      for(unsigned int j = 0; j < nhit; ++j){
        if(i == j) continue;

        const double dsq = util::sqr(x[i]-x[j])+util::sqr(y[i]-y[j]);
        closest = std::min(closest, dsq);
      }

      if(closest > util::sqr(fDHitGood*2)){
        w[i] = 0;
        ++ndrop;
      }
    }

    if(nhit-ndrop == 0) return 0;

    if(nhit-ndrop == 1){
      start[0] = x[0]-.1;
      end[0]   = x[0]+.1;
      start[1] = end[1] = y[0];
      return 1;
    }

    // check to see if there are only 2 points the fit is easy then
    if(nhit-ndrop == 2){
      start[0] = x[0];
      start[1] = y[0];
      end[0]   = x[1];
      end[1]   = y[1];
      if(start[0] > end[0]){
        std::swap(start[0], end[0]);
        std::swap(start[1], end[1]);
      }
      return 0;
    }

    bool droppedAny = true;
    while(droppedAny){
      chiSqr = geo::LinFitMinDperp(x, y, w,
                                   &start[0], &start[1],
                                   &end[0], &end[1]);

      LOG_DEBUG("CosmicTrack") << "end points from fit "
        << start[0] << " " << start[1] << " " << end[0] << " " << end[1];

      // loop over the hits and for the ones that are more than fDHitGood cm
      // off the line set the weight to be 0 so that they do not contribute to
      // future fits. Only drop the worst hit, and then iterate to get the next
      // worst and so on. Otherwise we can end up dropping lots of hits we
      // shouldn't due to one outlier.
      double dSqMax = 0;
      int dSqMaxIdx = -1;
      for(unsigned int n = 0; n < nhit; ++n){
        // Already dropped
        if(w[n] == 0) continue;

        const double dSq = geo::DsqrToLine(x[n], y[n],
                                           start[0], start[1],
                                           end[0], end[1]);
        if(dSq > dSqMax){
          dSqMax = dSq;
          dSqMaxIdx = n;
        }
      } // end for n

      droppedAny = false;
      if(dSqMax > util::sqr(fDHitGood)){
        LOG_DEBUG("CosmicTrack") << "dropping hit at ("
                                 << x[dSqMaxIdx] << ", " << y[dSqMaxIdx] << ")"
                                 << " dSq = " << dSqMax;

        w[dSqMaxIdx] = 0;
        droppedAny = true;
        ++ndrop;
        assert(nhit-ndrop >= 2); // If we're down to 2, we should be done
      }
    }// end while

    bool addedAny = false;
    // loop over the dropped points and see if any can be added back
    for(unsigned int n = 0; n < nhit; ++n){
      // Already part of the fit
      if(w[n] != 0) continue;

      const double dSq = geo::DsqrToLine(x[n], y[n],
                                         start[0], start[1],
                                         end[0], end[1]);
      if(dSq < util::sqr(fDHitGood)){
        LOG_DEBUG("CosmicTrack")
          << "adding back hit at (" << x[n] << ", " << y[n]
          << ")" << " dSq = " << dSq;

        w[n] = w_orig[n];
        --ndrop;
        addedAny = true;
      }
    } // end for n

    if(addedAny){
      // Fit with this final set of points
      chiSqr = geo::LinFitMinDperp(x, y, w,
                                   &start[0], &start[1],
                                   &end[0], &end[1]);
    }

    LOG_DEBUG("CosmicTrack") << "end points from fit "
                             << start[0] << " " << start[1] << " "
                             << end[0] << " " << end[1];

    // improve end point resolution by using distance of closest approach
    // ignore cases with completely horizontal/vertical lines to prevent failures/ infinities

    if(fabs(start[0]-end[0])>0.01&&fabs(start[1]-end[1])>0.01) {

      const TVector3 slope(end[0]-start[0],end[1]-start[1],0);
      const TVector3 intercept(0,(-end[0]/slope[0])*slope[1]+end[1],0); // define fit line
      double tmin[2], tmax[2];
      TVector3 closest; // will be point of closest approach on fit line
      tmin[0] = 9999.;
      tmax[0] = -9999.;
      tmin[1] = 9999.;
      tmax[1] = -9999.;

      for(unsigned int n = 0; n < nhit; ++n){
        if(w[n] == 0) continue;

        TVector3 point(x[n],y[n],0);

        geo::ClosestApproach(point,intercept,slope,closest);

        // save extremum; will be new end points
        if(closest[0]<tmin[0]) {  // since no vertical lines allowed, can just use extremum in z
          tmin[0] = closest[0];
          tmin[1] = closest[1];
        }
        if(closest[0]>tmax[0]) {
          tmax[0] = closest[0];
          tmax[1] = closest[1];
        }

        // check which end point each extremum is closer to, set new end points accordingly
        if(fabs(tmin[0]-end[0])<fabs(tmin[0]-start[0])) { 
          end[0] = tmin[0];
          end[1] = tmin[1];
          start[0] = tmax[0];
          start[1] = tmax[1];
        }
        else {
          end[0] = tmax[0];
          end[1] = tmax[1];
          start[0] = tmin[0];
          start[1] = tmin[1];
        }
      }
    }
    // Communicate back which points we dropped to the caller
    w_orig = w;

    return ndrop;
  }

  //......................................................................
  bool CosmicTrackAlg::IsTrackDownstreamFromTiming(const rb::Cluster* slice) const
  {
    art::ServiceHandle<geo::Geometry> geom;

    std::vector<double> zs, ts, ws;

    const int hitMax = slice->NCell();
    for(int hitIdx = 0; hitIdx < hitMax; ++hitIdx){
      art::Ptr<rb::CellHit> chit = slice->Cell(hitIdx);

      double xyz[3];
      geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
      zs.push_back(xyz[2]);
      ts.push_back(chit->TNS());
      ws.push_back(1);
    }

    // figure out the direction of the muon
    double ztStart[2] = {0.};
    double ztEnd[2]   = {0.};
    geo::LinFitMinDperp(zs, ts, ws,
                        &ztStart[0], &ztStart[1],
                        &ztEnd[0], &ztEnd[1]);

    // the slope should be positive if time increases with z, negative
    // otherwise
    const double ztSlope = (ztEnd[1]-ztStart[1])/(ztEnd[0]-ztStart[0]);

    return ztSlope > 0;
  }

}