Track3D.cxx

Go to the documentation of this file.

#include "HoughTrack/Track3D.h"

#include "Utilities/func/MathUtil.h"

#include "HoughTrack/Constants.h"
#include "Geometry/Geometry.h"
#include "RecoBaseHit/CellHit.h"

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

#include <TVector3.h>

#include <iostream>
#include <map>


htk::Track3D::Track3D()
{
}



htk::Track3D::Track3D(Track2D const& x_track, Track2D const& y_track) :
  default_value_(-9999),
  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_);
  bool yt_fit_valid = time_fit(y_track, dydt_, dzdt_y_track_);

  if (xt_fit_valid && yt_fit_valid)
  {
    // add a check to make sure it's almost the same for the two views
    //   having a large discrepancy in the two views might be a powerful
    //   discriminant

    dzdt_ = dzdt_x_track_;
    double velocity = util::pythag(dxdt_, dydt_, dzdt_);
    beta_ = velocity / Constants::SPEED_OF_LIGHT;
  } else {
    dzdt_ = default_value_;
    beta_ = default_value_;
  }
  
  for (auto const& x_hit : x_track.cluster.AllCells())
  {
      cluster_.Add(x_hit);
  }

  for (auto const& y_hit : y_track.cluster.AllCells())
  {
      cluster_.Add(y_hit);
  }
}



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



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



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



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



double htk::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 -9999;
}



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



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



bool htk::Track3D::time_fit
(Track2D const& track, double & dvdt, double & dzdt) 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}, {"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("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)
  {
  //  throw art::Exception(art::errors::InvalidNumber)
  //   << "Monopole/Track3D.cxx\n"
  //    << "Variance of t is zero.  Aborting.\n" << std::endl;
  
    dvdt = 999999;
    dzdt = 999999;
  }
  
  else
  {
    dvdt = variance.at("vt") / variance.at("t");
    dzdt = variance.at("zt") / variance.at("t");
  }

  //double v0 = mean.at("v") - dvdt * mean.at("t");
  //double z0 = mean.at("z") - dzdt * mean.at("t");
  
  return true;
}



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



double htk::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);
}