SimpleSegmentPairFilter.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/segmentPair/SimpleSegmentPairFilter.h>
 
#include <tracking/trackFindingCDC/fitting/CDCAxialStereoFusion.h>
 
#include <tracking/trackingUtilities/eventdata/tracks/CDCSegmentPair.h>
#include <tracking/trackingUtilities/eventdata/segments/CDCSegment2D.h>
#include <tracking/trackingUtilities/eventdata/trajectories/CDCTrajectory3D.h>
#include <tracking/trackingUtilities/eventdata/trajectories/CDCTrajectory2D.h>
 
#include <Math/VectorUtil.h>
#include <Math/Vector2D.h>
 
using namespace Belle2;
using namespace TrackFindingCDC;
using namespace TrackingUtilities;
 
Weight SimpleSegmentPairFilter::operator()(const CDCSegmentPair& segmentPair)
{
  const CDCSegment2D* ptrFromSegment = segmentPair.getFromSegment();
  const CDCSegment2D* ptrToSegment = segmentPair.getToSegment();
 
  assert(ptrFromSegment);
  assert(ptrToSegment);
 
  const CDCSegment2D& fromSegment = *ptrFromSegment;
  const CDCSegment2D& toSegment = *ptrToSegment;
 
  // Do fits
  const CDCTrajectory2D& fromFit = getFittedTrajectory2D(fromSegment);
  const CDCTrajectory2D& toFit = getFittedTrajectory2D(toSegment);
 
  // Check if segments are coaligned
  bool toSegmentIsCoaligned = fromFit.getTotalArcLength2D(toSegment) >= 0.0;
  bool fromSegmentIsCoaligned = toFit.getTotalArcLength2D(fromSegment) >= 0.0;
 
  if (not toSegmentIsCoaligned or not fromSegmentIsCoaligned) {
    return NAN;
  }
 
  // Check if there is a positive gap between from and to segment
  double fromFitGap = fromFit.getArcLength2DGap(fromSegment, toSegment);
  double toFitGap = toFit.getArcLength2DGap(fromSegment, toSegment);
 
  if (fromFitGap < -5 or fromFitGap > 50 or toFitGap < -5 or toFitGap > 50) {
    return NAN;
  }
 
  double fromFitFrontOffset = fromFit.getArcLength2DFrontOffset(fromSegment, toSegment);
  double toFitBackOffset = toFit.getArcLength2DBackOffset(fromSegment, toSegment);
 
  if (fromFitFrontOffset < 0 or
      fromFitFrontOffset > 50 or
      toFitBackOffset < 0 or
      toFitBackOffset > 50) {
    return NAN;
  }
 
  ROOT::Math::XYVector fromBackRecoPos2D = fromSegment.back().getRecoPos2D();
  ROOT::Math::XYVector toFrontRecoPos2D = toSegment.front().getRecoPos2D();
 
  // Momentum agreement cut
  ROOT::Math::XYVector fromMom2DAtFromBack = fromFit.getFlightDirection2D(fromBackRecoPos2D);
  ROOT::Math::XYVector toMom2DAtToFront = toFit.getFlightDirection2D(toFrontRecoPos2D);
 
  ROOT::Math::XYVector fromMom2DAtToFront = fromFit.getFlightDirection2D(toFrontRecoPos2D);
  ROOT::Math::XYVector toMom2DAtFromBack = toFit.getFlightDirection2D(fromBackRecoPos2D);
 
  double momAngleAtFromBack = ROOT::Math::VectorUtil::DeltaPhi(fromMom2DAtFromBack, toMom2DAtFromBack);
  double momAngleAtToFront = ROOT::Math::VectorUtil::DeltaPhi(toMom2DAtToFront, fromMom2DAtToFront);
 
  if (fabs(momAngleAtToFront) > 1.0 or fabs(momAngleAtFromBack) > 1.0) {
    return NAN;
  }
 
  // Proximity cut
  double fromFit_dist2DToFront_toSegment = fromFit.getDist2D(toSegment.front().getRecoPos2D());
  double toFit_dist2DToBack_fromSegment = toFit.getDist2D(fromSegment.back().getRecoPos2D());
 
  if (fromFit_dist2DToFront_toSegment < 10 and  toFit_dist2DToBack_fromSegment < 10) {
    getFittedTrajectory3D(segmentPair);
    return fromSegment.size() + toSegment.size();
  } else {
    return NAN;
  }
}
 
const CDCTrajectory2D& SimpleSegmentPairFilter::getFittedTrajectory2D(const CDCSegment2D& segment) const
{
  CDCTrajectory2D& trajectory2D = segment.getTrajectory2D();
  if (not trajectory2D.isFitted()) {
    getRiemannFitter().update(trajectory2D, segment);
  }
  return trajectory2D;
}
 
const CDCTrajectory3D& SimpleSegmentPairFilter::getFittedTrajectory3D(const CDCSegmentPair& segmentPair) const
{
  const CDCSegment2D* ptrFromSegment = segmentPair.getFromSegment();
  const CDCSegment2D* ptrToSegment = segmentPair.getToSegment();
 
  const CDCSegment2D& fromSegment = *ptrFromSegment;
  const CDCSegment2D& toSegment = *ptrToSegment;
 
  // Do fits if still necessary.
  getFittedTrajectory2D(fromSegment);
  getFittedTrajectory2D(toSegment);
 
  CDCAxialStereoFusion fusionFit;
  fusionFit.reconstructFuseTrajectories(segmentPair);
  return segmentPair.getTrajectory3D();
}