CDCTrajectory2D.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
#include <tracking/trackingUtilities/eventdata/trajectories/CDCTrajectory2D.h>
 
#include <tracking/trackingUtilities/eventdata/trajectories/CDCBFieldUtil.h>
 
#include <cdc/topology/CDCWireTopology.h>
#include <cdc/topology/CDCWireLayer.h>
#include <cdc/topology/WireLine.h>
#include <cdc/topology/ISuperLayer.h>
 
#include <tracking/trackingUtilities/geometry/UncertainPerigeeCircle.h>
#include <tracking/trackingUtilities/geometry/PerigeeCircle.h>
#include <tracking/trackingUtilities/geometry/VectorUtil.h>
 
#include <tracking/trackingUtilities/numerics/EForwardBackward.h>
#include <tracking/trackingUtilities/numerics/ESign.h>
#include <tracking/trackingUtilities/numerics/Quadratic.h>
 
#include <framework/gearbox/Const.h>
 
#include <vector>
#include <utility>
#include <ostream>
#include <cmath>
#include <cassert>
 
using namespace Belle2;
using namespace CDC;
using namespace TrackingUtilities;
 
CDCTrajectory2D::CDCTrajectory2D()
  : m_localOrigin()
  , m_localPerigeeCircle()
{
}
 
CDCTrajectory2D::CDCTrajectory2D(const UncertainPerigeeCircle& perigeeCircle)
  : m_localOrigin(0.0, 0.0)
  , m_localPerigeeCircle(perigeeCircle)
{
}
 
CDCTrajectory2D::CDCTrajectory2D(const ROOT::Math::XYVector& localOrigin,
                                 const UncertainPerigeeCircle& localPerigeeCircle,
                                 double flightTime)
  : m_localOrigin(localOrigin)
  , m_localPerigeeCircle(localPerigeeCircle)
  , m_flightTime(flightTime)
{
}
 
CDCTrajectory2D::CDCTrajectory2D(const ROOT::Math::XYVector& pos2D,
                                 const double time,
                                 const ROOT::Math::XYVector& mom2D,
                                 const double charge,
                                 const double bZ)
  : m_localOrigin(pos2D)
  , m_localPerigeeCircle(CDCBFieldUtil::absMom2DToCurvature(mom2D.R(), charge, bZ),
                         VectorUtil::unit(mom2D),
                         0.0)
  , m_flightTime(time)
{
}
 
CDCTrajectory2D::CDCTrajectory2D(const ROOT::Math::XYVector& pos2D,
                                 const double time,
                                 const ROOT::Math::XYVector& mom2D,
                                 const double charge)
  : m_localOrigin(pos2D)
  , m_localPerigeeCircle(CDCBFieldUtil::absMom2DToCurvature(mom2D.R(), charge, pos2D),
                         VectorUtil::unit(mom2D),
                         0.0)
  , m_flightTime(time)
{
}
 
bool CDCTrajectory2D::isFitted() const
{
  return not getLocalCircle()->isInvalid();
}
 
void CDCTrajectory2D::clear()
{
  m_localOrigin.SetXY(0.0, 0.0);
  m_localPerigeeCircle.invalidate();
  m_flightTime = NAN;
}
 
void CDCTrajectory2D::reverse()
{
  m_localPerigeeCircle.reverse();
  m_flightTime = -m_flightTime;
}
 
 
CDCTrajectory2D CDCTrajectory2D::reversed() const
{
  CDCTrajectory2D result = *this;
  result.reverse();
  return result;
}
 
std::array<double, 2> CDCTrajectory2D::reconstructBothZ(const WireLine& wireLine,
                                                        const double distance,
                                                        const double z) const
{
  ROOT::Math::XYVector globalPos2D = wireLine.sagPos2DAtZ(z);
  ROOT::Math::XYVector movePerZ = wireLine.sagMovePerZ(z);
 
  ROOT::Math::XYVector localPos2D = globalPos2D - getLocalOrigin();
  const PerigeeCircle& localCircle = getLocalCircle();
 
  double fastDistance = distance != 0.0 ? localCircle.fastDistance(distance) : 0.0;
 
  double c = localCircle.fastDistance(localPos2D) - fastDistance;
  double b = localCircle.gradient(localPos2D).Dot(movePerZ);
  double a = localCircle.n3() * movePerZ.Mag2();
 
  const std::pair<double, double> solutionsDeltaZ = solveQuadraticABC(a, b, c);
 
  // Put the solution of smaller deviation first
  const std::array<double, 2> solutionsZ{solutionsDeltaZ.second + z, solutionsDeltaZ.first + z};
  return solutionsZ;
}
 
double CDCTrajectory2D::reconstructZ(const WireLine& wireLine,
                                     const double distance,
                                     const double z) const
{
  const std::array<double, 2> solutionsZ = reconstructBothZ(wireLine, distance, z);
 
  bool firstIsInCDC = (wireLine.backwardZ() < solutionsZ[0] and
                       solutionsZ[0] < wireLine.forwardZ());
  bool secondIsInCDC = (wireLine.backwardZ() < solutionsZ[1] and
                        solutionsZ[1] < wireLine.forwardZ());
 
  // Prefer the solution with the smaller deviation from the given z position which is the first
  assert(not(std::fabs(solutionsZ[0] - z) > std::fabs(solutionsZ[1] - z)));
  const double recoZ = (firstIsInCDC or not secondIsInCDC) ? solutionsZ[0] : solutionsZ[1];
  return recoZ;
}
 
std::array<ROOT::Math::XYZVector, 2> CDCTrajectory2D::reconstructBoth3D(const WireLine& wireLine,
    const double distance,
    const double z) const
{
  const std::array<double, 2> solutionsZ = reconstructBothZ(wireLine, distance, z);
 
  const ROOT::Math::XYZVector firstRecoWirePos3D = wireLine.sagPos3DAtZ(solutionsZ[0]);
  const ROOT::Math::XYZVector secondRecoWirePos3D = wireLine.sagPos3DAtZ(solutionsZ[1]);
  const auto& tmp1 = getClosest(VectorUtil::getXYVector(firstRecoWirePos3D));
  const auto& tmp2 = getClosest(VectorUtil::getXYVector(secondRecoWirePos3D));
  return {{{tmp1.X(), tmp1.Y(), firstRecoWirePos3D.z()},
      {tmp2.X(), tmp2.Y(), secondRecoWirePos3D.z()}
    }};
}
 
ROOT::Math::XYZVector CDCTrajectory2D::reconstruct3D(const WireLine& wireLine,
                                                     const double distance,
                                                     const double z) const
{
  const double recoZ = reconstructZ(wireLine, distance, z);
  const ROOT::Math::XYZVector recoWirePos2D = wireLine.sagPos3DAtZ(recoZ);
  const auto& tmp = getClosest(VectorUtil::getXYVector(recoWirePos2D));
  return ROOT::Math::XYZVector(tmp.X(), tmp.Y(), recoZ);
}
 
ROOT::Math::XYVector CDCTrajectory2D::getClosest(const ROOT::Math::XYVector& point) const
{
  return getLocalCircle()->closest(point - getLocalOrigin()) + getLocalOrigin();
}
 
ISuperLayer CDCTrajectory2D::getISuperLayerAfter(ISuperLayer iSuperLayer, bool movingOutward) const
{
  if (ISuperLayerUtil::isInvalid(iSuperLayer)) return ISuperLayerUtil::c_Invalid;
 
  ISuperLayer minimalISuperLayer = getMinimalISuperLayer();
  ISuperLayer maximalISuperLayer = getMaximalISuperLayer();
  if (minimalISuperLayer == maximalISuperLayer) return ISuperLayerUtil::c_Invalid; // No next super layer to go to
  if (iSuperLayer == minimalISuperLayer) return ISuperLayerUtil::getNextOutwards(iSuperLayer);
  if (iSuperLayer == maximalISuperLayer) return ISuperLayerUtil::getNextInwards(iSuperLayer);
 
  if (movingOutward) {
    return ISuperLayerUtil::getNextOutwards(iSuperLayer);
  } else {
    return ISuperLayerUtil::getNextInwards(iSuperLayer);
  }
}
 
ISuperLayer CDCTrajectory2D::getISuperLayerAfterStart(bool movingOutward) const
{
  ISuperLayer iSuperLayer = getStartISuperLayer();
  return getISuperLayerAfter(iSuperLayer, movingOutward);
}
 
ISuperLayer CDCTrajectory2D::getISuperLayerAfterStart(const EForwardBackward forwardBackwardInfo) const
{
  bool movingOutward = isMovingOutward();
  if (forwardBackwardInfo == EForwardBackward::c_Backward) {
    movingOutward = not movingOutward;
  }
  return getISuperLayerAfterStart(movingOutward);
}
 
ISuperLayer CDCTrajectory2D::getNextISuperLayer() const
{
  return getISuperLayerAfterStart(EForwardBackward::c_Forward);
}
 
ISuperLayer CDCTrajectory2D::getPreviousISuperLayer() const
{
  return getISuperLayerAfterStart(EForwardBackward::c_Backward);
}
 
ISuperLayer CDCTrajectory2D::getAxialISuperLayerAfterStart(const EForwardBackward forwardBackwardInfo) const
{
  bool movingOutward = isMovingOutward();
  if (forwardBackwardInfo == EForwardBackward::c_Backward) {
    movingOutward = not movingOutward;
  }
  ISuperLayer startISuperLayer = getStartISuperLayer();
  if (ISuperLayerUtil::isInvalid(startISuperLayer)) return ISuperLayerUtil::c_Invalid;
 
  ISuperLayer nextISuperLayer = getISuperLayerAfter(startISuperLayer, movingOutward);
  if (ISuperLayerUtil::isInvalid(nextISuperLayer)) return ISuperLayerUtil::c_Invalid;
  if (ISuperLayerUtil::isAxial(nextISuperLayer)) return nextISuperLayer;
 
  ISuperLayer iSuperLayerStep = nextISuperLayer - startISuperLayer;
  assert(std::abs(iSuperLayerStep) == 1);
  bool nextMovingOutward = iSuperLayerStep > 0;
  return getISuperLayerAfter(nextISuperLayer, nextMovingOutward);
}
 
ISuperLayer CDCTrajectory2D::getNextAxialISuperLayer() const
{
  return getAxialISuperLayerAfterStart(EForwardBackward::c_Forward);
}
 
ISuperLayer CDCTrajectory2D::getPreviousAxialISuperLayer() const
{
  return getAxialISuperLayerAfterStart(EForwardBackward::c_Backward);
}
 
ISuperLayer CDCTrajectory2D::getMaximalISuperLayer() const
{
  double maximalCylindricalR = getMaximalCylindricalR();
  return CDCWireTopology::getInstance().getISuperLayerAtCylindricalR(maximalCylindricalR);
}
 
ISuperLayer CDCTrajectory2D::getStartISuperLayer() const
{
  double startCylindricalR = getLocalOrigin().R();
  return CDCWireTopology::getInstance().getISuperLayerAtCylindricalR(startCylindricalR);
}
 
ISuperLayer CDCTrajectory2D::getMinimalISuperLayer() const
{
  double minimalCylindricalR = getMinimalCylindricalR();
  return CDCWireTopology::getInstance().getISuperLayerAtCylindricalR(minimalCylindricalR);
}
 
bool CDCTrajectory2D::isCurler(double factor) const
{
  const CDCWireTopology& topology = CDCWireTopology::getInstance();
  return getMaximalCylindricalR() < factor * topology.getOuterCylindricalR();
}
 
bool CDCTrajectory2D::isOriginer(double factor) const
{
  const CDCWireTopology& topology = CDCWireTopology::getInstance();
  return getMinimalCylindricalR() < factor * topology.getInnerCylindricalR();
}
 
 
ESign CDCTrajectory2D::getChargeSign() const
{
  return CDCBFieldUtil::ccwInfoToChargeSign(getLocalCircle()->orientation());
}
 
double CDCTrajectory2D::getAbsMom2D(const double bZ) const
{
  return CDCBFieldUtil::curvatureToAbsMom2D(getLocalCircle()->curvature(), bZ);
}
 
double CDCTrajectory2D::getAbsMom2D() const
{
  ROOT::Math::XYVector position = getSupport();
  return CDCBFieldUtil::curvatureToAbsMom2D(getLocalCircle()->curvature(), position);
}
 
ROOT::Math::XYVector CDCTrajectory2D::getInnerExit() const
{
  const CDCWireTopology& topology = CDCWireTopology::getInstance();
  const CDCWireLayer& innerMostLayer = topology.getWireLayers().front();
  double innerCylindricalR = innerMostLayer.getInnerCylindricalR();
 
  const ROOT::Math::XYVector support = getSupport();
  const PerigeeCircle globalCircle = getGlobalCircle();
  if (support.Mag2() < innerCylindricalR * innerCylindricalR) {
    // If we start within the inner volume of the CDC we want the trajectory to enter the CDC
    // and not stop at first intersection with the inner wall.
    // Therefore we take the inner exit that comes after the apogee (far point of the circle).
    const ROOT::Math::XYVector apogee = globalCircle.apogee();
    return globalCircle.atCylindricalRForwardOf(apogee, innerCylindricalR);
 
  } else {
    return globalCircle.atCylindricalRForwardOf(support, innerCylindricalR);
  }
}
 
ROOT::Math::XYVector CDCTrajectory2D::getOuterExit(double factor) const
{
  const CDCWireTopology& topology = CDCWireTopology::getInstance();
  const CDCWireLayer& outerMostLayer = topology.getWireLayers().back();
  double outerCylindricalR = outerMostLayer.getOuterCylindricalR() * factor;
 
  const ROOT::Math::XYVector support = getSupport();
  const PerigeeCircle globalCircle = getGlobalCircle();
  if (support.Mag2() > outerCylindricalR * outerCylindricalR) {
    // If we start outside of the volume of the CDC we want the trajectory to enter the CDC
    // and not stop at first intersection with the outer wall.
    // Therefore we take the outer exit that comes after the perigee.
    const ROOT::Math::XYVector perigee = globalCircle.perigee();
    return globalCircle.atCylindricalRForwardOf(perigee, outerCylindricalR);
 
  } else {
    return getGlobalCircle().atCylindricalRForwardOf(support, outerCylindricalR);
  }
}
 
ROOT::Math::XYVector CDCTrajectory2D::getExit() const
{
  const ROOT::Math::XYVector outerExit = getOuterExit();
  const ROOT::Math::XYVector innerExit = getInnerExit();
  const ROOT::Math::XYVector localExit =  getLocalCircle()->chooseNextForwardOf(ROOT::Math::XYVector(0, 0),
                                          outerExit - getLocalOrigin(),
                                          innerExit - getLocalOrigin());
  return localExit + getLocalOrigin();
}
 
void CDCTrajectory2D::setPosMom2D(const ROOT::Math::XYVector& pos2D,
                                  const ROOT::Math::XYVector& mom2D,
                                  const double charge)
{
  m_localOrigin = pos2D;
  double curvature = CDCBFieldUtil::absMom2DToCurvature(mom2D.R(), charge, pos2D);
  ROOT::Math::XYVector phiVec = VectorUtil::unit(mom2D);
  double impact = 0.0;
  m_localPerigeeCircle = UncertainPerigeeCircle(curvature, phiVec, impact);
}
 
double CDCTrajectory2D::setLocalOrigin(const ROOT::Math::XYVector& localOrigin)
{
  double arcLength2D = calcArcLength2D(localOrigin);
  m_flightTime += arcLength2D / Const::speedOfLight;
  m_localPerigeeCircle.passiveMoveBy(localOrigin - m_localOrigin);
  m_localOrigin = localOrigin;
  return arcLength2D;
}
 
 
std::ostream& TrackingUtilities::operator<<(std::ostream& output, const CDCTrajectory2D& trajectory2D)
{
  return output << "Local origin : " << trajectory2D.getLocalOrigin() << ", "
         << "local circle : " << trajectory2D.getLocalCircle();
}