Prong.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    Prong.cxx
// \brief   A Prong is a Cluster with defined start position and direction
// \version $Id: Prong.cxx,v 1.12 2012-09-21 02:32:54 greenc Exp $
// \author  $Author: greenc $
// \date    $Date: 2012-09-21 02:32:54 $
////////////////////////////////////////////////////////////////////////
#include "RecoBase/Prong.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/CellHit.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"

namespace rb
{
  //......................................................................
  // 3D
  Prong::Prong(const art::PtrVector<rb::CellHit>& cells,
         TVector3 start, TVector3 dir, int id)
    : rb::Cluster(cells, id)
  {
    SetStart(start);
    SetDir(dir);
  }

  Prong::Prong(const rb::Cluster& clust,
         TVector3 start, TVector3 dir, int id)
    : rb::Cluster(clust)
  {
    SetID(id);
    SetStart(start);
    SetDir(dir);
  }

  //......................................................................
  // 2D prong
  Prong::Prong(const art::PtrVector<rb::CellHit>& cells,
         geo::View_t view, double v0, double z0, double dv, double dz, int id)
    : rb::Cluster(view, cells, id)
  {
    SetStart(v0, z0);
    SetDir(dv, dz);
  }

  //......................................................................
  // 2D prong
  Prong::Prong(const std::vector< art::Ptr<rb::CellHit> >& cells,
               geo::View_t view, double v0, double z0, double dv, double dz, int id)
    : rb::Cluster(view, cells, id)
  {
    SetStart(v0, z0);
    SetDir(dv, dz);
  }

  Prong::Prong(const rb::Cluster& clust,
         double v0, double z0, double dv, double dz, int id)
    : rb::Cluster(clust)
  {
    SetID(id);
    SetStart(v0, z0);
    SetDir(dv, dz);
  }

  //......................................................................

  Prong::~Prong()
  {
  }

  //......................................................................

  void Prong::SetStart(TVector3 start)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    // I've removed this assert to allow for 2D prongs to set a 3D start
    // point if they so choose.
    // assert(fView == geo::kXorY);

    fStart = start;
  }

  //......................................................................

  void Prong::SetStart(double v0, double z0)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    assert(fView == geo::kX || fView == geo::kY);

    if(fView == geo::kX){
      fStart = TVector3(v0, 0, z0);
    }
    else{
      fStart = TVector3(0, v0, z0);
    }
  }

  //......................................................................

  void Prong::SetDir(TVector3 dir)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    // I've removed this assert to allow for 2D prongs to set a 3D start
    // direction if they so choose.
    // assert(fView == geo::kXorY);

    fDir = dir.Unit(); // Ensure normalized
  }

  //......................................................................

  void Prong::SetDir(double dv, double dz)
  {
    fPrecalcTotalGeV = -1; // Invalidate any cached value

    assert(fView == geo::kX || fView == geo::kY);

    if(fView == geo::kX){
      fDir = TVector3(dv, 0, dz);
    }
    else{
      fDir = TVector3(0, dv, dz);
    }

    fDir = fDir.Unit(); // Ensure normalized
  }

  //......................................................................

  double Prong::W(const rb::CellHit* chit) const
  {
   
    //TODO: If the prong is 2D, return a position of 0
    //for the opposite view.  This is in line with what happens
    //in the w method for a track and cluster.  In the future
    //trying to do this should fail in all cases.  This patch
    //keeps anything from breaking until there is a long-term fix.

    // Removing this line so that 2D prongs can get a sensible W value.
    // if (chit->View() == fView) return 0;

    art::ServiceHandle<geo::Geometry> geom;

    double xyz[3];
    geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
    const double z1 = xyz[2];

    /// \todo CJB - we're just moving the point in the XY plane to get
    /// it on the line. Maybe we should take the route perpendicular
    /// to the line instead. Similarly for Track

    const double dt = (z1-fStart.Z())/fDir.Z();

    if(chit->View() == geo::kX) return fStart.Y()+dt*fDir.Y();
    if(geom->DetId() == novadaq::cnv::kTESTBEAM) return -(fStart.X()+dt*fDir.X()); // TB horizontal modules (Y-view) have readout on opposite end as ND/FD
    return fStart.X()+dt*fDir.X();

    /*
    // Solve the problem of the closest approach of the trajectory to the line
    // made by the cell by projecting down to a 2D trajectory and 2D point for
    // the cell.
    TVector3 start2D = fStart;
    TVector3 dir2D = fDir;
    TVector3 pt2D(xyz);
    const int otherView = 1-chit->View();
    start2D[otherView] = dir2D[otherView] = pt2D[otherView] = 0;

    // Actually find the closest point
    TVector3 closest;
    geo::ClosestApproach(pt2D, start2D, dir2D, closest);

    // Need to solve for the W position of the line, the exact part that we
    // projected away above. So just figure out Z and evaluate the prong there.
    const double dt = (closest.Z()-fStart.Z())/fDir.Z();

    return (fStart+dt*fDir)[otherView];
    */
  }

  //......................................................................
  double Prong::TotalLength() const
  {
    art::ServiceHandle<geo::Geometry> geom;

    double ret = 0;

    const unsigned int cellMax = NCell();

    for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
      const art::Ptr<rb::CellHit> chit = Cell(cellIdx);

      TVector3 xyz;
      geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);

      const TVector3 dxyz = (chit->View() == geo::kX) ? TVector3(0, 1, 0)
                                                      : TVector3(1, 0, 0);

      double dist;
      geo::ClosestApproach(fStart, fStart+fDir,
                           xyz, xyz+dxyz,
                           &dist);

      ret = std::max(ret, dist);
    }
    return ret;
  }

  //......................................................................
  double Prong::DistanceFromStart(double z) const
  {
    // A prong is a straight line so this is easy

    // For a completely vertical track any answer is nonsense
    if(fDir.Z() == 0) return 0;

    return (z-fStart.Z())/fDir.Z();
  }

  //......................................................................

  bool Prong::operator< (const Prong& other) const
  {
    return this->Start().Z() < other.Start().Z();
  }

} // end namespace rb
////////////////////////////////////////////////////////////////////////