VXDMomentumEstimationTools.h

/**************************************************************************
 * 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.                  *
 **************************************************************************/
 
#pragma once
 
#include <tracking/spacePointCreation/SpacePoint.h>
#include <vxd/dataobjects/VxdID.h>
#include <framework/geometry/B2Vector3.h>
#include <framework/dataobjects/Helix.h>
#include <vxd/geometry/GeoCache.h>
#include <vxd/geometry/SensorInfoBase.h>
#include <algorithm>
#include <cmath>
#include <framework/logging/Logger.h>
 
 
namespace Belle2 {
  class SVDCluster;
  class PXDCluster;

  template <class ClusterType>
  class VXDMomentumEstimationTools {
 
  private:
    VXDMomentumEstimationTools() { }
    VXDMomentumEstimationTools(const VXDMomentumEstimationTools&);
    VXDMomentumEstimationTools& operator = (const VXDMomentumEstimationTools&);
 
  public:
    static const VXDMomentumEstimationTools& getInstance()
    {
      static VXDMomentumEstimationTools instance;
      return instance;
    }

    double getDEDX(const ClusterType& cluster, const ROOT::Math::XYZVector& momentum, const ROOT::Math::XYZVector& position,
                   short charge) const
    {
 
      const Helix trajectory(position, momentum, charge, 1.5);
      const double calibratedCharge = getCalibratedCharge(cluster);
      const double pathLength = getPathLength(cluster, trajectory);
 
      return calibratedCharge / pathLength;
    }

    static double getDEDXWithThickness(const ClusterType& cluster)
    {
      const double calibratedCharge = getCalibratedCharge(cluster);
      const double pathLength = getThicknessOfCluster(cluster);
 
      return calibratedCharge / pathLength;
    }

    static double getThicknessOfCluster(const ClusterType& cluster)
    {
      const VxdID&   vxdID = cluster.getSensorID();
      const VXD::SensorInfoBase& sensorInfoBase = VXD::GeoCache::getInstance().getSensorInfo(vxdID);
      return sensorInfoBase.getThickness();
    }

    static double getWidthOfCluster(const ClusterType& cluster)
    {
      const VxdID&   vxdID = cluster.getSensorID();
      const VXD::SensorInfoBase& sensorInfoBase = VXD::GeoCache::getInstance().getSensorInfo(vxdID);
      return sensorInfoBase.getWidth();
    }

    static double getLengthOfCluster(const ClusterType& cluster)
    {
      const VxdID&   vxdID = cluster.getSensorID();
      const VXD::SensorInfoBase& sensorInfoBase = VXD::GeoCache::getInstance().getSensorInfo(vxdID);
      return sensorInfoBase.getLength();
    }

    double getRadiusOfCluster(const ClusterType& cluster) const
    {
      const SpacePoint* spacePoint = cluster.template getRelated<SpacePoint>("SpacePoints");
 
      double radius = 0;
      if (spacePoint == nullptr) {
        // Fallback function
        VxdID vxdID = cluster.getSensorID();
        VxdID::baseType layer = vxdID.getLayerNumber();
        radius = m_layerPositions[layer - 1] / 100.0;
        B2WARNING("Could not determine SpacePoint for this cluster. Falling back to geometrical information.");
      } else {
        // The positions is the one with w = 0 which is the one on the lowest detector plane.
        radius = spacePoint->getPosition().Perp();
      }
      return radius;
    }

    static VxdID::baseType getLayerOfCluster(const ClusterType& cluster)
    {
      VxdID vxdID = cluster.getSensorID();
      VxdID::baseType layer = vxdID.getLayerNumber();
      return layer;
    }

    static VxdID::baseType getLadderOfCluster(const ClusterType& cluster)
    {
      VxdID vxdID = cluster.getSensorID();
      VxdID::baseType ladder = vxdID.getLadderNumber();
      return ladder;
    }

    static VxdID::baseType getSensorNumberOfCluster(const ClusterType& cluster)
    {
      VxdID vxdID = cluster.getSensorID();
      VxdID::baseType sensorNumber = vxdID.getSensorNumber();
      return sensorNumber;
    }

    static VxdID::baseType getSegmentNumberOfCluster(const ClusterType& cluster)
    {
      VxdID vxdID = cluster.getSensorID();
      VxdID::baseType segmentNumber = vxdID.getSegmentNumber();
      return segmentNumber;
    }

    static double getCalibratedCharge(const ClusterType& cluster)
    {
      const double charge = cluster.getCharge();
      const double calibration = getCalibration();
 
      return calibration * charge;
    }

    ROOT::Math::XYZVector getEntryMomentumOfMCParticle(const ClusterType&) const
    {
      B2FATAL("Can not deal with this cluster type!");
      return ROOT::Math::XYZVector();
    }

    ROOT::Math::XYZVector getEntryPositionOfMCParticle(const ClusterType&) const
    {
      B2FATAL("Can not deal with this cluster type!");
      return ROOT::Math::XYZVector();
    }

    double getPathLength(const ClusterType& cluster, const Helix& trajectory) const
    {
      // This is not quite correct but we can not do better without doing an extrapolation with material effects.
      // If the distance_3D is nan, it means that the helix does not reach into the cluster or curls that much
      // that it does not reach the far end of the cluster. The first case is strange
      // because the track is associated with the cluster (so it should normally reach it). The second case should be
      // really rare (because the clusters are that thin), but we have to deal with it.
      // There is also the possibility that the track curls that much that it interacts with the same sensor twice.
      // This can not be handled properly in this stage.
 
      const double thickness = getThicknessOfCluster(cluster);
      const double radius = getRadiusOfCluster(cluster);
      const double width = getWidthOfCluster(cluster);
      const double length = getLengthOfCluster(cluster);
 
      const double perp_s_at_cluster_entry = trajectory.getArcLength2DAtCylindricalR(radius);
 
      if (std::isnan(perp_s_at_cluster_entry)) {
        // This case is really strange. I do not know how to deal with it probably.
        return thickness;
      }
 
      const ROOT::Math::XYZVector& position_at_inner_radius = trajectory.getPositionAtArcLength2D(perp_s_at_cluster_entry);
 
      const double perp_s_at_cluster_exit = trajectory.getArcLength2DAtCylindricalR(radius + thickness);
 
      if (std::isnan(perp_s_at_cluster_exit)) {
        return std::min(width, length);
      }
 
      const ROOT::Math::XYZVector& position_at_outer_radius = trajectory.getPositionAtArcLength2D(perp_s_at_cluster_exit);
 
      const double distance_3D = (position_at_outer_radius - position_at_inner_radius).R();
 
      return distance_3D;
    }
 
  private:
    const double m_layerPositions[6] = {1.42, 2.18, 3.81, 8.0, 10.51, 13.51};

    static double getCalibration()
    {
      return 1;
    }
  };

  template<>
  double VXDMomentumEstimationTools<PXDCluster>::getCalibration();

  template <>
  ROOT::Math::XYZVector VXDMomentumEstimationTools<PXDCluster>::getEntryMomentumOfMCParticle(const PXDCluster& cluster) const;

  template <>
  ROOT::Math::XYZVector VXDMomentumEstimationTools<SVDCluster>::getEntryMomentumOfMCParticle(const SVDCluster& cluster) const;

  template <>
  ROOT::Math::XYZVector VXDMomentumEstimationTools<PXDCluster>::getEntryPositionOfMCParticle(const PXDCluster& cluster) const;

  template <>
  ROOT::Math::XYZVector VXDMomentumEstimationTools<SVDCluster>::getEntryPositionOfMCParticle(const SVDCluster& cluster) const;
}