TrackExtrap.h

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.

/// \file TrackExtrap.h
/// \brief Definition of tools for track extrapolation
///
/// \author Philippe Pillot, Subatech

#ifndef O2_MCH_TRACKEXTRAP_H_
#define O2_MCH_TRACKEXTRAP_H_

#include <cstddef>
#include <optional>

#include <TMatrixD.h>

namespace o2
{
namespace mch
{

class TrackParam;

/// Class holding tools for track extrapolation
class TrackExtrap
{
 public:
  // static class
  TrackExtrap() = delete;
  ~TrackExtrap() = delete;

  TrackExtrap(const TrackExtrap&) = delete;
  TrackExtrap& operator=(const TrackExtrap&) = delete;
  TrackExtrap(TrackExtrap&&) = delete;
  TrackExtrap& operator=(TrackExtrap&&) = delete;

  static void setField();

  /// Return true if the field is switched ON
  static bool isFieldON() { return sFieldON; }

  /// Switch to Runge-Kutta extrapolation v2
  static void useExtrapV2(bool extrapV2 = true) { sExtrapV2 = extrapV2; }

  static double getImpactParamFromBendingMomentum(double bendingMomentum);
  static double getBendingMomentumFromImpactParam(double impactParam);

  static void linearExtrapToZ(TrackParam& trackParam, double zEnd);
  static void linearExtrapToZCov(TrackParam& trackParam, double zEnd, bool updatePropagator = false);

  static bool extrapToZ(TrackParam& trackParam, double zEnd);
  static bool extrapToZCov(TrackParam& trackParam, double zEnd, bool updatePropagator = false);

  static bool extrapToVertex(TrackParam& trackParam, double xVtx, double yVtx, double zVtx, double errXVtx, double errYVtx)
  {
    /// Extrapolate track parameters to vertex, corrected for multiple scattering and energy loss effects
    /// Add branson correction resolution and energy loss fluctuation to parameter covariances
    return extrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, true, true);
  }
  static bool extrapToVertexWithoutELoss(TrackParam& trackParam, double xVtx, double yVtx, double zVtx, double errXVtx, double errYVtx)
  {
    /// Extrapolate track parameters to vertex, corrected for multiple scattering effects only
    /// Add branson correction resolution to parameter covariances
    return extrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, true, false);
  }
  static bool extrapToVertexWithoutBranson(TrackParam& trackParam, double zVtx,
                                           double xUpstream = 0., double yUpstream = 0.,
                                           std::optional<double> zUpstream = std::nullopt)
  {
    /// Extrapolate track parameters to vertex, corrected for energy loss effects only
    /// Add dispersion due to multiple scattering and energy loss fluctuation to parameter covariances
    return extrapToVertex(trackParam, 0., 0., zVtx, 0., 0., false, true, xUpstream, yUpstream, zUpstream);
  }
  static bool extrapToVertexUncorrected(TrackParam& trackParam, double zVtx,
                                        double xUpstream = 0., double yUpstream = 0.,
                                        std::optional<double> zUpstream = std::nullopt)
  {
    /// Extrapolate track parameters to vertex without multiple scattering and energy loss corrections
    /// Add dispersion due to multiple scattering to parameter covariances
    return extrapToVertex(trackParam, 0., 0., zVtx, 0., 0., false, false, xUpstream, yUpstream, zUpstream);
  }

  static bool extrapToMID(TrackParam& trackParam);

  static double getMCSAngle2(const TrackParam& param, double dZ, double x0);
  static void addMCSEffect(TrackParam& trackParam, double dZ, double x0);

  static void printNCalls();

 private:
  static bool extrapToVertex(TrackParam& trackParam, double xVtx, double yVtx, double zVtx,
                             double errXVtx, double errYVtx, bool correctForMCS, bool correctForEnergyLoss,
                             double xUpstream = 0., double yUpstream = 0., std::optional<double> zUpstream = std::nullopt);

  static bool getAbsorberCorrectionParam(double trackXYZIn[3], double trackXYZOut[3], double pTotal,
                                         double& pathLength, double& f0, double& f1, double& f2,
                                         double& meanRho, double& totalELoss, double& sigmaELoss2);

  static void addMCSEffectInAbsorber(TrackParam& param, double signedPathLength, double f0, double f1, double f2);

  static double betheBloch(double pTotal, double pathLength, double rho, double atomicZ, double atomicZoverA);
  static double energyLossFluctuation(double pTotal, double pathLength, double rho, double atomicZoverA);

  static bool correctMCSEffectInAbsorber(TrackParam& param, double xVtx, double yVtx, double zVtx, double errXVtx, double errYVtx,
                                         double absZBeg, double pathLength, double f0, double f1, double f2);
  static void correctELossEffectInAbsorber(TrackParam& param, double eLoss, double sigmaELoss2);

  static void cov2CovP(const TMatrixD& param, TMatrixD& cov);
  static void covP2Cov(const TMatrixD& param, TMatrixD& covP);

  static void convertTrackParamForExtrap(TrackParam& trackParam, double forwardBackward, double* v3);
  static void recoverTrackParam(double* v3, double Charge, TrackParam& trackParam);

  static bool extrapToZRungekutta(TrackParam& trackParam, double zEnd);
  static bool extrapToZRungekuttaV2(TrackParam& trackParam, double zEnd);
  static bool extrapOneStepRungekutta(double charge, double step, const double* vect, double* vout);

  static constexpr double SMuMass = 0.105658;                         ///< Muon mass (GeV/c2)
  static constexpr double SAbsZBeg = -90.;                            ///< Position of the begining of the absorber (cm)
  static constexpr double SAbsZEnd = -505.;                           ///< Position of the end of the absorber (cm)
  static constexpr double SSimpleBPosition = -0.5 * (994.05 + 986.6); ///< Position of the dipole (cm)
  static constexpr double SSimpleBLength = 0.5 * (502.1 + 309.4);     ///< Length of the dipole (cm)
  static constexpr int SMaxStepNumber = 5000;                         ///< Maximum number of steps for track extrapolation
  static constexpr double SRungeKuttaMaxResidue = 0.002;              ///< Max z-distance to destination to stop the track extrap (cm)
  static constexpr double SRungeKuttaMaxResidueV2 = 0.01;             ///< Max z-distance to destination to stop the track extrap v2 (cm)
  /// Most probable value (GeV/c) of muon momentum in bending plane (used when B = 0)
  /// Needed to get some "reasonable" corrections for MCS and E loss even if B = 0
  static constexpr double SMostProbBendingMomentum = 2.;
  static constexpr double SMuonFilterZBeg = -1471.;    ///< Position of the begining of the muon filter (cm)
  static constexpr double SMuonFilterThickness = 120.; ///< Thickness of the muon filter (cm)
  /// Position of the end of the muon filter (cm)
  static constexpr double SMuonFilterZEnd = SMuonFilterZBeg - SMuonFilterThickness;
  static constexpr double SMuonFilterX0 = 1.76; ///< Radiation length of the muon filter (cm)
  static constexpr double SMIDZ = -1603.5;      ///< Position of the first MID chamber (cm)

  static bool sExtrapV2; ///< switch to Runge-Kutta extrapolation v2

  static double sSimpleBValue; ///< Magnetic field value at the centre
  static bool sFieldON;        ///< true if the field is switched ON

  static std::size_t sNCallExtrapToZCov; ///< number of times the method extrapToZCov(...) is called
  static std::size_t sNCallField;        ///< number of times the method Field(...) is called
};

} // namespace mch
} // namespace o2

#endif // O2_MCH_TRACKEXTRAP_H_