SimplePXDStateFilter.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/ckf/pxd/filters/states/SimplePXDStateFilter.h>
 
#include <tracking/trackingUtilities/eventdata/trajectories/CDCTrajectory3D.h>
#include <tracking/trackingUtilities/eventdata/trajectories/CDCTrajectory2D.h>
#include <tracking/trackingUtilities/eventdata/trajectories/CDCTrajectorySZ.h>
 
#include <tracking/trackingUtilities/eventdata/trajectories/CDCBFieldUtil.h>
#include <tracking/spacePointCreation/SpacePoint.h>
#include <tracking/dataobjects/RecoTrack.h>
 
#include <Math/Vector2D.h>
#include <Math/Vector3D.h>
 
using namespace Belle2;
using namespace TrackingUtilities;
 
namespace {
  unsigned int getPTRange(const ROOT::Math::XYZVector& momentum)
  {
    const double pT = momentum.Rho();
    if (pT > 0.4) {
      return 0;
    } else if (pT > 0.2) {
      return 1;
    } else {
      return 2;
    }
  }
}
 
constexpr const SimplePXDStateFilter::MaximalValueArray SimplePXDStateFilter::m_param_maximumHelixDistanceXY;
constexpr const SimplePXDStateFilter::MaximalValueArray SimplePXDStateFilter::m_param_maximumResidual;
constexpr const SimplePXDStateFilter::MaximalValueArray SimplePXDStateFilter::m_param_maximumChi2;
 
void SimplePXDStateFilter::beginRun()
{
  m_cachedBField = TrackingUtilities::CDCBFieldUtil::getBFieldZ();
}
 
Weight SimplePXDStateFilter::operator()(const BasePXDStateFilter::Object& pair)
{
  const std::vector<TrackingUtilities::WithWeight<const CKFToPXDState*>>& previousStates = pair.first;
  CKFToPXDState* currentState = pair.second;
 
  const auto* spacePoint = currentState->getHit();
 
  genfit::MeasuredStateOnPlane firstMeasurement;
  if (currentState->mSoPSet()) {
    firstMeasurement = currentState->getMeasuredStateOnPlane();
  } else {
    firstMeasurement = previousStates.back()->getMeasuredStateOnPlane();
  }
 
  ROOT::Math::XYZVector position = ROOT::Math::XYZVector(firstMeasurement.getPos());
  ROOT::Math::XYZVector momentum = ROOT::Math::XYZVector(firstMeasurement.getMom());
 
  const ROOT::Math::XYZVector hitPosition = spacePoint->getPosition();
 
  const bool sameHemisphere = fabs(position.Phi() - hitPosition.Phi()) < TMath::PiOver2();
  if (not sameHemisphere) {
    return NAN;
  }
 
  double valueToCheck;
  const MaximalValueArray* maximumValues;
 
  if (not currentState->isFitted() and not currentState->mSoPSet()) {
    // Filter 1
    const RecoTrack* cdcTrack = previousStates.front()->getSeed();
 
    B2ASSERT("A path without a seed?", cdcTrack);
 
    const CDCTrajectory3D trajectory(position, 0, momentum, cdcTrack->getChargeSeed(), m_cachedBField);
 
    const double arcLength = trajectory.calcArcLength2D(hitPosition);
    const ROOT::Math::XYVector& trackPositionAtHit2D = trajectory.getTrajectory2D().getPos2DAtArcLength2D(arcLength);
    const double trackPositionAtHitZ = trajectory.getTrajectorySZ().mapSToZ(arcLength);
    const ROOT::Math::XYZVector trackPositionAtHit(trackPositionAtHit2D.X(), trackPositionAtHit2D.Y(), trackPositionAtHitZ);
    const ROOT::Math::XYZVector differenceHelix = trackPositionAtHit - hitPosition;
 
    valueToCheck = differenceHelix.Rho();
    maximumValues = &m_param_maximumHelixDistanceXY;
  } else if (not currentState->isFitted()) {
    // Filter 2
    const double residual = m_kalmanStepper.calculateResidual(firstMeasurement, *currentState);
 
    valueToCheck = residual;
    maximumValues = &m_param_maximumResidual;
  } else {
    // Filter 3
    valueToCheck = currentState->getChi2();
    maximumValues = &m_param_maximumChi2;
  }
 
  const unsigned int layer = currentState->getGeometricalLayer();
  if (valueToCheck > (*maximumValues)[layer - 1][getPTRange(momentum)]) {
    return NAN;
  }
 
  return valueToCheck;
}