Track3D.cxx

Go to the documentation of this file.

#include "Monopole/Track3D.h"

#include "Monopole/Constants.h"
#include "Geometry/Geometry.h"

#include <art/Framework/Services/Registry/ServiceHandle.h>
// #include <canvas/Utilities/Exception.h>
#include <art/Utilities/Exception.h>

#include <TVector3.h>

#include <algorithm>
#include <iostream>
#include <map>



mono::Track3D::Track3D()
{
}



mono::Track3D::Track3D(Track2D const& x_track, Track2D const& y_track) :
  default_value_  (-9999),
  dxdt_           (default_value_),
  dydt_           (default_value_),
  dzdt_           (default_value_),
  dzdt_x_track_   (default_value_),
  dzdt_y_track_   (default_value_),
  r2_xt_          (default_value_),
  r2_yt_          (default_value_),
  r2_zt_x_track_  (default_value_),
  r2_zt_y_track_  (default_value_),
  max_time_gap_xz_(default_value_),
  max_time_gap_yz_(default_value_),
  beta_           (default_value_),
  x_track_        (x_track),
  y_track_        (y_track)
{
  x0_ = x_track.intercept;
  y0_ = y_track.intercept;
  z0_ = 0.0;

  dxdz_ = x_track.slope;
  dydz_ = y_track.slope;

  bool xt_fit_valid =
    time_fit(x_track, dxdt_, dzdt_x_track_, r2_xt_, r2_zt_x_track_);
  bool yt_fit_valid =
    time_fit(y_track, dydt_, dzdt_y_track_, r2_yt_, r2_zt_y_track_);

  if (xt_fit_valid && yt_fit_valid)
  {
    dzdt_ = dzdt_x_track_;
    double velocity = util::pythag(dxdt_, dydt_, dzdt_);
    beta_ = velocity / Constants::SPEED_OF_LIGHT;
  }


  /*
    Hit Matching

    Note that this can be done here or at the beginning of the constructor,
    but it does not affect the x_track and y_track themselves, which are used
    for the fit.
  */
  for (auto const& x_hit : x_track.cluster.AllCells())
  {
    bool has_match = false;
    for (auto const& y_hit : y_track.cluster.AllCells())
      if (abs(int(x_hit->Plane()) - int(y_hit->Plane())) <= 1)
        has_match = true;
    
    if (has_match)
      cluster_.Add(x_hit);
  }

  for (auto const& y_hit : y_track.cluster.AllCells())
  {
    bool has_match = false;
    for (auto const& x_hit : x_track.cluster.AllCells())
      if (abs(int(x_hit->Plane()) - int(y_hit->Plane())) <= 1)
        has_match = true;

    if (has_match)
      cluster_.Add(y_hit);
  }

  max_time_gap_xz_ = time_gap(cluster_.XCells());
  max_time_gap_yz_ = time_gap(cluster_.YCells());
}



TVector3 mono::Track3D::start() const
{
  return TVector3(x0_, y0_, z0_);
}



TVector3 mono::Track3D::velocity() const
{
  return TVector3(dxdt_, dydt_, dzdt_);
}



double mono::Track3D::dxdz() const
{
  return dxdz_;
}



double mono::Track3D::dydz() const
{
  return dydz_;
}



double mono::Track3D::dzdt(geo::View_t view) const
{
  if (view == geo::View_t::kX)
    return dzdt_x_track_;

  if (view == geo::View_t::kY)
    return dzdt_y_track_;

  return default_value_;
}



double mono::Track3D::r2_xt() const
{
  return r2_xt_;
}



double mono::Track3D::r2_yt() const
{
  return r2_yt_;
}



double mono::Track3D::r2_zt(geo::View_t view) const
{
  if (view == geo::View_t::kX)
    return r2_zt_x_track_;

  if (view == geo::View_t::kY)
    return r2_zt_y_track_;

  return default_value_;
}



double mono::Track3D::max_time_gap_xz() const
{
  return max_time_gap_xz_;
}



double mono::Track3D::max_time_gap_yz() const
{
  return max_time_gap_yz_;
}



double mono::Track3D::beta() const
{
  return beta_;
}



rb::Cluster mono::Track3D::cluster() const
{
  return cluster_;
}



bool mono::Track3D::time_fit
(Track2D const& track,
 double & dvdt, double & dzdt,
 double & r2_vt, double & r2_zt) const
{
  dvdt = default_value_;
  dzdt = default_value_;
  double N = track.cluster.NCell();
  // we need at least two hits to get a non-zero variance
  if (N <= 2)
    return false;
  
  art::ServiceHandle<geo::Geometry> geometry;

  std::map<std::string, double> mean = { {"v", 0.0}, {"z", 0.0}, {"t", 0.0} };
  for (auto const& hit : track.cluster.AllCells())
  {
    double xyz[3];
    geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
    TVector3 hit_vector(xyz);
    
    mean.at("v") += geometry->CellTpos(hit->Plane(), hit->Cell()) / N;
    mean.at("z") += hit_vector.z() / N;
    mean.at("t") += hit->TNS() / N;
  }

  std::map<std::string, double> variance =
    { {"t", 0.0}, {"v", 0.0}, {"z", 0.0}, {"vt", 0.0}, {"zt", 0.0} };
  for (auto const& hit : track.cluster.AllCells())
  {
    double xyz[3];
    geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
    TVector3 hit_vector(xyz);
    double v = geometry->CellTpos(hit->Plane(), hit->Cell());
    double z = hit_vector.z();
    double t = hit->TNS();

    variance.at("t")  += get_variance(t, mean.at("t")) / (N - 1);
    variance.at("v")  += get_variance(v, mean.at("v")) / (N - 1);
    variance.at("z")  += get_variance(z, mean.at("z")) / (N - 1);
    variance.at("vt") +=
      get_covariance(v, mean.at("v"), t, mean.at("t")) / (N - 1);
    variance.at("zt") +=
      get_covariance(z, mean.at("z"), t, mean.at("t")) / (N - 1);
  }
  
  if (variance.at("t") == 0  || variance.at("v") == 0 || variance.at("z") == 0)
    return false;
  // throw art::Exception(art::errors::InvalidNumber)
  //   << "Monopole/Track3D.cxx\n"
  //   << "Variance is zero.  Aborting.\n" << std::endl;


  //
  // Set all of the results:
  //
  dvdt = variance.at("vt") / variance.at("t");
  dzdt = variance.at("zt") / variance.at("t");

  // The intercept is currently not used.
  // double v0 = mean.at("v") - dvdt * mean.at("t");
  // double z0 = mean.at("z") - dzdt * mean.at("t");

  // Calculate the squared correlation coefficienct (r^2):
  r2_vt =
    std::pow(variance.at("vt"), 2) / (variance.at("v") * variance.at("t"));
  r2_zt =
    std::pow(variance.at("zt"), 2) / (variance.at("z") * variance.at("t"));

  return true;
}



double mono::Track3D::time_gap(art::PtrVector<rb::CellHit> const& hits) const
{
  double result = default_value_;

  std::set<double> times;
  for (auto const& hit : hits)
    times.insert(hit->TNS());

  std::adjacent_find
  (times.begin(), times.end(), [&](double object, double next_object)
  {
    double gap = next_object - object;
    result = std::max(gap, result);

    // This return is meaningless, but a boolean return is required by
    // adjacent find.  We choose false, so the loop goes until the end.
    return false;
  });

  return result;
}



double mono::Track3D::get_variance
(double const& a, double const& mean_a) const
{
  return (a - mean_a) * (a - mean_a);
}



double mono::Track3D::get_covariance
(double const& a, double const& mean_a,
 double const& b, double const& mean_b) const
{
  return (a - mean_a) * (b - mean_b);
}



std::ostream& mono::operator<<(std::ostream& os, mono::Track3D const& t)
{
  art::ServiceHandle<geo::Geometry> geometry;

  os << "Track3D: dxdz, dydz = " << t.dxdz() << ", " << t.dydz() << "\n";

  for (auto const& hit : t.cluster().AllCells())
  {
    std::string view;
    if (hit->View() == geo::View_t::kX)
      view = "XZ";
    else if (hit->View() == geo::View_t::kY)
      view = "YZ";
    else
      view = "BAD";

    double xyz[3];
    geometry->CellInfo(hit->Plane(), hit->Cell(), NULL, xyz);
    TVector3 hit_vector(xyz);
    double v = geometry->CellTpos(hit->Plane(), hit->Cell());
    double z = hit_vector.z();
    double t = hit->TNS();
    
    os << "Track3D: View " << view
       << " Hit: t, v, z  = " << t << ", " << v << ", " << z << "\n";
  }
  
  return os;
}