InteractPlane.cxx

Go to the documentation of this file.

#include <cmath>
#include "lardata/RecoObjects/InteractPlane.h"
#include "lardata/RecoObjects/SurfPlane.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "cetlib_except/exception.h"

namespace trkf {

  InteractPlane::InteractPlane(double tcut) :
    Interactor(tcut)
  {}

  InteractPlane::~InteractPlane()
  {}

  bool InteractPlane::noise(const KTrack& trk, double s, TrackError& noise_matrix) const
  {
    // Get LAr service.

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

    // Make sure we are on a plane surface (throw exception if not).

    const SurfPlane* psurf = dynamic_cast<const SurfPlane*>(&*trk.getSurface());
    if(psurf == 0)
      throw cet::exception("InteractPlane") 
        << "InteractPlane called for non-planar surface.\n";

    // Clear noise matrix.

    noise_matrix.clear();

    // Unpack track parameters.

    const TrackVector& vec = trk.getVector();
    double dudw = vec[2];
    double dvdw = vec[3];
    double pinv = vec[4];
    double mass = trk.Mass();

    // If distance is zero, or momentum is infinite, return zero noise.

    if(pinv == 0. || s == 0.)
      return true;

    // Make a crude estimate of the range of the track.

    double p = 1./std::abs(pinv);
    double p2 = p*p;
    double e2 = p2 + mass*mass;
    double e = std::sqrt(e2);
    double t = e - mass;
    double dedx = 0.001 * detprop->Eloss(p, mass, getTcut());
    double range = t / dedx;
    if(range > 100.)
      range = 100.;

    // Calculate the radiation length in cm.

    double x0 = larprop->RadiationLength() / detprop->Density();

    // Calculate projected rms scattering angle.
    // Use the estimted range in the logarithm factor.
    // Use the incremental propagation distance in the square root factor.

    double betainv = std::sqrt(1. + pinv*pinv * mass*mass);
    double theta_fact = (0.0136 * pinv * betainv) * (1. + 0.038 * std::log(range/x0));
    double theta02 = theta_fact*theta_fact * std::abs(s/x0);

    // Calculate some sommon factors needed for multiple scattering.

    double ufact2 = 1. + dudw*dudw;
    double vfact2 = 1. + dvdw*dvdw;
    double uvfact2 = 1. + dudw*dudw + dvdw*dvdw;
    double uvfact = std::sqrt(uvfact2);
    double uv = dudw * dvdw;
    double dist2_3 = s*s / 3.;
    double dist_2 = std::abs(s) / 2.;
    if(trk.getDirection() == Surface::BACKWARD)
      dist_2 = -dist_2;

    // Calculate energy loss fluctuations.

    double evar = 1.e-6 * detprop->ElossVar(p, mass) * std::abs(s); // E variance (GeV^2).
    double pinvvar = evar * e2 / (p2*p2*p2);                        // Inv. p variance (1/GeV^2)

    // Fill elements of noise matrix.

    // Position submatrix.

    noise_matrix(0,0) = dist2_3 * theta02 * ufact2;           // sigma^2(u,u)
    noise_matrix(1,0) = dist2_3 * theta02 * uv;               // sigma^2(u,v)
    noise_matrix(1,1) = dist2_3 * theta02 * vfact2;           // sigma^2(v,v)

    // Slope submatrix.

    noise_matrix(2,2) = theta02 * uvfact2 * ufact2;           // sigma^2(u', u')
    noise_matrix(3,2) = theta02 * uvfact2 * uv;               // sigma^2(v', u')
    noise_matrix(3,3) = theta02 * uvfact2 * vfact2;           // sigma^2(v', v')

    // Same-view position-slope correlations.

    noise_matrix(2,0) = dist_2 * theta02 * uvfact * ufact2;   // sigma^2(u', u)
    noise_matrix(3,1) = dist_2 * theta02 * uvfact * vfact2;   // sigma^2(v', v)

    // Opposite-view position-slope correlations.

    noise_matrix(2,1) = dist_2 * theta02 * uvfact * uv;       // sigma^2(u', v)
    noise_matrix(3,0) = dist_2 * theta02 * uvfact * uv;       // sigma^2(v', u)

    // Momentum correlations (zero).

    noise_matrix(4,0) = 0.;                                   // sigma^2(pinv, u)
    noise_matrix(4,1) = 0.;                                   // sigma^2(pinv, v)
    noise_matrix(4,2) = 0.;                                   // sigma^2(pinv, u')
    noise_matrix(4,3) = 0.;                                   // sigma^2(pinv, v')

    // Energy loss fluctuations.

    noise_matrix(4,4) = pinvvar;                              // sigma^2(pinv, pinv)

    // Done (success).

    return true;
  }

} // end namespace trkf