Propagator.cxx

Go to the documentation of this file.

#include "lardata/RecoObjects/Propagator.h"
#include "lardata/RecoObjects/SurfXYZPlane.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "cetlib_except/exception.h"

namespace trkf {

  Propagator::Propagator(double tcut, bool doDedx,
                         const std::shared_ptr<const Interactor>& interactor) :
    fTcut(tcut),
    fDoDedx(doDedx),
    fInteractor(interactor)
  {}

  Propagator::~Propagator()
  {}

  boost::optional<double> Propagator::vec_prop(KTrack& trk,
                                               const std::shared_ptr<const Surface>& psurf, 
                                               PropDirection dir,
                                               bool doDedx,
                                               TrackMatrix* prop_matrix,
                                               TrackError* noise_matrix) const
  {
    // Default result.
    
    auto result = boost::make_optional<double>(false, 0.);

    // Get the inverse momentum (assumed to be track parameter four).

    double pinv = trk.getVector()(4);

    // If dE/dx is not requested, or if inverse momentum is zero, then
    // it is safe to propagate in one step.  In this case, just pass
    // the call to short_vec_prop with unlimited distance.

    bool dedx = getDoDedx() && doDedx;
    if(!dedx || pinv == 0.)
      result = short_vec_prop(trk, psurf, dir, dedx, prop_matrix, noise_matrix);

    else {

      // dE/dx is requested.  In this case we limit the maximum
      // propagation distance such that the kinetic energy of the
      // particle should not change by more thatn 10%.

      // Get LAr service.

      auto const * detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();

      // Initialize propagation matrix to unit matrix (if specified).

      int nvec = trk.getVector().size();
      if(prop_matrix)
        *prop_matrix = ublas::identity_matrix<TrackVector::value_type>(nvec);

      // Initialize noise matrix to zero matrix (if specified).

      if(noise_matrix) {
        noise_matrix->resize(nvec, nvec, false);
        noise_matrix->clear();
      }

      // Remember the starting track.

      KTrack trk0(trk);

      // Make pointer variables pointing to local versions of the
      // propagation and noise matrices, or null if not specified.

      TrackMatrix local_prop_matrix;
      TrackMatrix* plocal_prop_matrix = (prop_matrix==0 ? 0 : &local_prop_matrix);
      TrackError local_noise_matrix;
      TrackError* plocal_noise_matrix = (noise_matrix==0 ? 0 : &local_noise_matrix);

      // Cumulative propagation distance.

      double s = 0.;
        
      // Begin stepping loop.
      // We put a maximum iteration count to prevent infinite loops caused by
      // floating point pathologies.  The iteration count is large enough to reach
      // any point in the tpc using the minimum step size (for a reasonable tpc).

      bool done = false;
      int nitmax = 10000;  // Maximum number of iterations.
      int nit = 0;         // Iteration count.
      while(!done) {

        // If the iteration count exceeds the maximum, return failure.

        ++nit;
        if(nit > nitmax) {
          trk = trk0;
          result = boost::optional<double>(false, 0.);
          return result;
        }

        // Estimate maximum step distance according to the above
        // stated principle.

        pinv = trk.getVector()(4);
        double mass = trk.Mass();
        double p = 1./std::abs(pinv);
        double e = std::hypot(p, mass);
        double t = p*p / (e + mass);
        double dedx = 0.001 * detprop->Eloss(p, mass, fTcut);
        double smax = 0.1 * t / dedx;
        if (smax <= 0.)
          throw cet::exception("Propagator") << __func__ << ": maximum step " << smax << "\n";

        // Always allow a step of at least 0.3 cm (about one wire spacing).

        if(smax < 0.3)
          smax = 0.3;

        // First do a test propagation (without dE/dx and errors) to
        // find the distance to the destination surface.

        KTrack trktest(trk);
        boost::optional<double> dist = short_vec_prop(trktest, psurf, dir, false, 0, 0);

        // If the test propagation failed, return failure.

        if(!dist) {
          trk = trk0;
          return dist;
        }

        // Generate destionation surface for this step (either final
        // destination, or some intermediate surface).

        std::shared_ptr<const Surface> pstep;
        if(std::abs(*dist) <= smax) {
          done = true;
          pstep = psurf;
        }
        else {

          // Generate intermediate surface.
          // First get point where track will intersect intermediate surface.

          double xyz0[3];                                // Starting point.
          trk.getPosition(xyz0);
          double xyz1[3];                                // Destination point.
          trktest.getPosition(xyz1);
          double frac = smax / std::abs(*dist);
          double xyz[3];                                 // Intermediate point.
          xyz[0] = xyz0[0] + frac * (xyz1[0] - xyz0[0]);
          xyz[1] = xyz0[1] + frac * (xyz1[1] - xyz0[1]);
          xyz[2] = xyz0[2] + frac * (xyz1[2] - xyz0[2]);

          // Choose orientation of intermediate surface perpendicular
          // to track.

          double mom[3];
          trk.getMomentum(mom);

          // Make intermediate surface object.

          pstep = std::shared_ptr<const Surface>(new SurfXYZPlane(xyz[0], xyz[1], xyz[2],
                                                                  mom[0], mom[1], mom[2]));
        }

        // Do the actual step propagation.

        dist = short_vec_prop(trk, pstep, dir, doDedx, 
                              plocal_prop_matrix, plocal_noise_matrix);

        // If the step propagation failed, return failure.

        if(!dist) {
          trk = trk0;
          return dist;
        }

        // Update cumulative propagation distance.

        s += *dist;

        // Update cumulative propagation matrix (left-multiply).

        if(prop_matrix != 0) {
          TrackMatrix temp = prod(*plocal_prop_matrix, *prop_matrix);
          *prop_matrix = temp;
        }

        // Update cumulative noise matrix.

        if(noise_matrix != 0) {
          TrackMatrix temp = prod(*noise_matrix, trans(*plocal_prop_matrix));
          TrackMatrix temp2 = prod(*plocal_prop_matrix, temp);
          *noise_matrix = ublas::symmetric_adaptor<TrackMatrix>(temp2);
          *noise_matrix += *plocal_noise_matrix;
        }
      }

      // Set the final result (distance + success).

      result = boost::optional<double>(true, s);
    }

    // Done.

    return result;
  }

  boost::optional<double> Propagator::lin_prop(KTrack& trk,
                                               const std::shared_ptr<const Surface>& psurf, 
                                               PropDirection dir,
                                               bool doDedx,
                                               KTrack* ref,
                                               TrackMatrix* prop_matrix,
                                               TrackError* noise_matrix) const
  {
    // Default result.

    boost::optional<double> result;

    if(ref == 0)
      result = vec_prop(trk, psurf, dir, doDedx, prop_matrix, noise_matrix);
    else {

      // A reference track has been provided.

      // It is an error (throw exception) if the reference track and
      // the track to be propagted are not on the same surface.

      if(!trk.getSurface()->isEqual(*(ref->getSurface())))
        throw cet::exception("Propagator") << 
          "Input track and reference track not on same surface.\n";

      // Remember the starting track and reference track.

      KTrack trk0(trk);
      KTrack ref0(*ref);

      // Propagate the reference track.  Make sure we calculate the
      // propagation matrix.

      TrackMatrix prop_temp;
      if(prop_matrix == 0)
        prop_matrix = &prop_temp;

      // Do the propgation.  The returned result will be the result of
      // this propagatrion.

      result = vec_prop(*ref, psurf, dir, doDedx, prop_matrix, noise_matrix);
      if(!!result) {

        // Propagation of reference track succeeded.  Update the track
        // state vector and surface of the track to be propagated.

        TrackVector diff = trk.getSurface()->getDiff(trk.getVector(), ref0.getVector());
        TrackVector newvec = ref->getVector() + prod(*prop_matrix, diff);

        // Store updated state vector and surface.

        trk.setVector(newvec);
        trk.setSurface(psurf);
        trk.setDirection(ref->getDirection());
        if (!trk.getSurface()->isEqual(*(ref->getSurface())))
          throw cet::exception("Propagator") << __func__ << ": surface mismatch";

        // Final validity check.  In case of failure, restore the track
        // and reference track to their starting values.

        if(!trk.isValid()) {
          result = boost::optional<double>(false, 0.);
          trk = trk0;
          *ref = ref0;
        }
      }
      else {

        // Propagation failed.
        // Restore the reference track to its starting value, so that we ensure 
        // the reference track and the actual track remain on the same surface.

        trk = trk0;
        *ref = ref0;
      }
    }

    // Done.

    return result;
  }


  boost::optional<double> Propagator::err_prop(KETrack& tre,
                                               const std::shared_ptr<const Surface>& psurf, 
                                               PropDirection dir,
                                               bool doDedx,
                                               KTrack* ref,
                                               TrackMatrix* prop_matrix) const
  {
    // Propagate without error, get propagation matrix.

    TrackMatrix prop_temp;
    if(prop_matrix == 0)
      prop_matrix = &prop_temp;
    boost::optional<double> result = lin_prop(tre, psurf, dir, doDedx, ref, prop_matrix, 0);

    // If propagation succeeded, update track error matrix.

    if(!!result) {
      TrackMatrix temp = prod(tre.getError(), trans(*prop_matrix));
      TrackMatrix temp2 = prod(*prop_matrix, temp);
      TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
      tre.setError(newerr);
    }

    // Done.

    return result;
  }

  boost::optional<double> Propagator::noise_prop(KETrack& tre,
                                                 const std::shared_ptr<const Surface>& psurf, 
                                                 PropDirection dir,
                                                 bool doDedx,
                                                 KTrack* ref) const
  {
    // Propagate without error, get propagation matrix and noise matrix.

    TrackMatrix prop_matrix;
    TrackError noise_matrix;
    boost::optional<double> result = lin_prop(tre, psurf, dir, doDedx, ref,
                                              &prop_matrix, &noise_matrix);

    // If propagation succeeded, update track error matrix.

    if(!!result) {
      TrackMatrix temp = prod(tre.getError(), trans(prop_matrix));
      TrackMatrix temp2 = prod(prop_matrix, temp);
      TrackError newerr = ublas::symmetric_adaptor<TrackMatrix>(temp2);
      newerr += noise_matrix;
      tre.setError(newerr);
    }

    // Done.

    return result;
  }

  boost::optional<double> Propagator::dedx_prop(double pinv, double mass,
                                                double s, double* deriv) const
  {
    // For infinite initial momentum, return with success status,
    // still infinite momentum.

    if(pinv == 0.)
      return boost::optional<double>(true, 0.);

    // Set the default return value to be uninitialized with value 0.

    boost::optional<double> result(false, 0.);

    // Get LAr service.

    auto const * detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();

    // Calculate final energy.

    double p1 = 1./std::abs(pinv);
    double e1 = std::hypot(p1, mass);
    double de = -0.001 * s * detprop->Eloss(p1, mass, fTcut);
    double emid = e1 + 0.5 * de;
    if(emid > mass) {
      double pmid = std::sqrt(emid*emid - mass*mass);
      double e2 = e1 - 0.001 * s * detprop->Eloss(pmid, mass, fTcut);
      if(e2 > mass) {
        double p2 = std::sqrt(e2*e2 - mass*mass);
        double pinv2 = 1./p2;
        if(pinv < 0.)
          pinv2 = -pinv2;

        // Calculation was successful, update result.

        result = boost::optional<double>(true, pinv2);

        // Also calculate derivative, if requested.

        if(deriv != 0)
          *deriv = pinv2*pinv2*pinv2 * e2 / (pinv*pinv*pinv * e1);
      }
    }

    // Done.

    return result;
  }
} // end namespace trkf