RtChebyshev.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
#include "MdtCalibData/RtChebyshev.h"
#include "MuonCalibMath/ChebychevPoly.h"
#include "GeoModelKernel/throwExcept.h"
using namespace MuonCalib;
 
RtChebyshev::RtChebyshev(const ParVec& vec) : 
    IRtRelation(vec) { 
    // check for consistency //
    if (nPar() < 3) {
        THROW_EXCEPTION("RtChebyshev::_init() - Not enough parameters!");
    }
    if (tLower() >= tUpper()) {
        THROW_EXCEPTION("Lower time boundary >= upper time boundary!");
    }
}  // end RtChebyshev::_init
 
std::string RtChebyshev::name() const { return "RtChebyshev"; }
double RtChebyshev::tBinWidth() const { return s_tBinWidth; }
 
double RtChebyshev::radius(double t) const {
    // INITIAL TIME CHECK //
    if (t < tLower()) return 0.0;
    if (t > tUpper()) return 14.6;
 
    // VARIABLES //
    // argument of the Chebyshev polynomials
    double x = getReducedTime(t);
    double rad{0.0};  // auxiliary radius
 
    // CALCULATE r(t) //
    for (unsigned int k = 0; k < nPar() - 2; k++) { 
        rad += par(k+2) * chebyshevPoly1st(k, x); 
    }
    return std::max(rad, 0.);
}
 
//*****************************************************************************
double RtChebyshev::driftVelocity(double t) const { 
    return (radius(t + 1.0) - radius(t));
    // Set derivative to 0 outside of the bounds
    if (t < tLower() || t > tUpper()) return 0.0;
 
    // Argument of the Chebyshev polynomials
    const double x = getReducedTime(t);
    // Chain rule
    const double dx_dt = 2. / (tUpper() - tLower());
    double drdt{0.};
    for (unsigned int k = 1; k < nPar() - 2; ++k) {
        // Calculate the contribution to dr/dt using k * U_{k-1}(x) * dx/dt
        drdt += par(k+2) *  chebyshevPoly1stPrime(k, x) * dx_dt;
    }
    return drdt; 
}
double RtChebyshev::driftAcceleration(double t) const {
    double acc{0.};
    // Argument of the Chebyshev polynomials
    const double x = getReducedTime(t);
    const double dx_dt = std::pow(2. / (tUpper() - tLower()), 2);
    for (unsigned int k = 2; k < nPar() - 2; ++k) {
        acc += par(k+2) *  chebyshevPoly1st2Prime(k, x) * dx_dt;
    }
    return acc * t;
}
double RtChebyshev::tLower() const { return par(0); }
double RtChebyshev::tUpper() const { return par(1); }
unsigned int RtChebyshev::numberOfRtParameters() const { return nPar() - 2; }
 
std::vector<double> RtChebyshev::rtParameters() const {
    return std::vector<double>{parameters().begin() +2, parameters().end()};
}
double RtChebyshev::getReducedTime(const double  t) const {
    return 2. * (t - 0.5 * (tUpper() + tLower())) / (tUpper() - tLower());
}