SimpleAxialSegmentPairFilter.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/trackFindingCDC/filters/axialSegmentPair/SimpleAxialSegmentPairFilter.h>
 
#include <tracking/trackingUtilities/eventdata/tracks/CDCAxialSegmentPair.h>
#include <tracking/trackingUtilities/eventdata/segments/CDCSegment2D.h>
 
#include <Math/Vector2D.h>
#include <Math/VectorUtil.h>
 
using namespace Belle2;
using namespace TrackFindingCDC;
using namespace TrackingUtilities;
 
SimpleAxialSegmentPairFilter::SimpleAxialSegmentPairFilter() : m_riemannFitter()
{
  m_riemannFitter.useOnlyOrientation();
}
 
Weight SimpleAxialSegmentPairFilter::operator()(const CDCAxialSegmentPair& axialSegmentPair)
{
  const CDCAxialSegment2D* ptrStartSegment = axialSegmentPair.getStartSegment();
  const CDCAxialSegment2D* ptrEndSegment = axialSegmentPair.getEndSegment();
 
  assert(ptrStartSegment);
  assert(ptrEndSegment);
 
  const CDCAxialSegment2D& startSegment = *ptrStartSegment;
  const CDCAxialSegment2D& endSegment = *ptrEndSegment;
 
  // Do fits
  const CDCTrajectory2D& startFit = getFittedTrajectory2D(startSegment);
  const CDCTrajectory2D& endFit = getFittedTrajectory2D(endSegment);
 
  // Check if segments are coaligned
  bool endSegmentIsCoaligned = startFit.getTotalArcLength2D(endSegment) >= 0.0;
  bool startSegmentIsCoaligned = endFit.getTotalArcLength2D(startSegment) >= 0.0;
 
  if (not endSegmentIsCoaligned or not startSegmentIsCoaligned) {
    return NAN;
  }
 
  // Check if there is a positive gap between start and end segment
  double startFitGap = startFit.getArcLength2DGap(startSegment, endSegment);
  double endFitGap = endFit.getArcLength2DGap(startSegment, endSegment);
 
  if (startFitGap < 0 or startFitGap > 100 or endFitGap < 0 or endFitGap > 100) {
    return NAN;
  }
 
  double startFitFrontOffset = startFit.getArcLength2DFrontOffset(startSegment, endSegment);
  double endFitBackOffset = endFit.getArcLength2DBackOffset(startSegment, endSegment);
 
  if (startFitFrontOffset < 0 or endFitBackOffset < 0) {
    return NAN;
  }
 
  ROOT::Math::XYVector startBackRecoPos2D = startSegment.back().getRecoPos2D();
  ROOT::Math::XYVector endFrontRecoPos2D = endSegment.front().getRecoPos2D();
 
  // Momentum agreement cut
  ROOT::Math::XYVector startMom2DAtStartBack = startFit.getFlightDirection2D(startBackRecoPos2D);
  ROOT::Math::XYVector endMom2DAtEndFront = endFit.getFlightDirection2D(endFrontRecoPos2D);
 
  ROOT::Math::XYVector startMom2DAtEndFront = startFit.getFlightDirection2D(endFrontRecoPos2D);
  ROOT::Math::XYVector endMom2DAtStartBack = endFit.getFlightDirection2D(startBackRecoPos2D);
 
  double momAngleAtStartBack = ROOT::Math::VectorUtil::DeltaPhi(startMom2DAtStartBack, endMom2DAtStartBack);
  double momAngleAtEndFront = ROOT::Math::VectorUtil::DeltaPhi(endMom2DAtEndFront, startMom2DAtEndFront);
 
  if (fabs(momAngleAtEndFront) > 2.0 or fabs(momAngleAtStartBack) > 2.0) {
    return NAN;
  }
 
  // Proximity cut
  double startFit_dist2DToFront_endSegment = startFit.getDist2D(endFrontRecoPos2D);
  double endFit_dist2DToBack_startSegment = endFit.getDist2D(startBackRecoPos2D);
 
  if (startFit_dist2DToFront_endSegment < 6 and endFit_dist2DToBack_startSegment < 6) {
    return startSegment.size() + endSegment.size();
  } else {
    return NAN;
  }
 
}
 
 
 
const CDCTrajectory2D& SimpleAxialSegmentPairFilter::getFittedTrajectory2D(const CDCAxialSegment2D& segment) const
{
 
  CDCTrajectory2D& trajectory2D = segment.getTrajectory2D();
  if (not trajectory2D.isFitted()) {
    getRiemannFitter().update(trajectory2D, segment);
  }
  return trajectory2D;
 
}
 
 
 
const CDCTrajectory2D& SimpleAxialSegmentPairFilter::getFittedTrajectory2D(const CDCAxialSegmentPair&
    axialSegmentPair) const
{
  CDCTrajectory2D& trajectory2D = axialSegmentPair.getTrajectory2D();
  if (not trajectory2D.isFitted()) {
    getRiemannFitter().update(trajectory2D, axialSegmentPair);
  }
  return trajectory2D;
}