MuonCalib::BFieldCorFunc Class Reference

#include <BFieldCorFunc.h>

Inheritance diagram for MuonCalib::BFieldCorFunc:

Collaboration diagram for MuonCalib::BFieldCorFunc:

Public Types
using	ParVec = std::vector< double >

Public Member Functions
	BFieldCorFunc (const std::string &quality, const CalibFunc::ParVec &parameters, const IRtRelation *rt)
	Constructor: quality = "high", slow but accurate initialization initialization of the correction function, quality = "medium", compromise between speed and accuracy of the initialization of the correction function (default), quality = "low", fast initialization of the correction function at the price of lower quality. More...
	BFieldCorFunc (const CalibFunc::ParVec &parameters, const IRtRelation *rt)
	Constructor: parameters[0] = high voltage [V], parameters[1] = $\epsilon$,. More...
double	epsilon () const
	< get the $\epsilon$ parameter of the B-field correction function More...
void	setRtRelationship (const IRtRelation &rt)
	< set the $\epsilon$ parameter of the B-field correction function = eps More...
std::string	name () const
	get the class name More...
double	correction (double t, double B_wire, double B_mu) const
	get t(r, $\vec{B}$!=0)-t(r, $\vec{B}$=0); t = drift time t [ns] for B=0; B_wire = magnetic field parallel to the anode wire of the given tube, B_mu = magnetic field orthogonal to the wire and parallel to the muon trajectory in the given tube [B] = Tesla More...
double	correction_to_B (double t, double B_wire, double B_mu, double B_factor=-1.0) const
virtual std::string	typeName () const
unsigned int	nPar () const
const ParVec &	parameters () const
double	par (unsigned int index) const

Private Member Functions
void	init (const std::string &quality, const CalibFunc::ParVec &params, const IRtRelation *rt)
double	t_from_r (const double r, const IRtRelation *rt) const
double	integral (const double r_min, const double r_max, const IRtRelation *rt) const

Private Attributes
std::vector< double >	m_param {}
std::string	m_quality {}
double	m_step_size {0.}
Amg::VectorX	m_alpha
double	m_t_min {-std::numeric_limits<double>::max()}
double	m_t_max {std::numeric_limits<double>::max()}
double	m_r_min {-std::numeric_limits<double>::max()}
double	m_r_max {std::numeric_limits<double>::max()}
ParVec	m_parameters {}

Detailed Description

This class allows the user to get the difference between the drift time measured by a tube operated in a magnetic field $\vec{B}$ and the drift time which would be measured by this tube if $\vec{B}$ vanished.

Correction:

$$t(r,\vec{B}) = t(r,\vec{B}=0) + B_\perp^{2-\epsilon}\cdot \int\limits_{25\ \mu m}^{r} \frac{v_{B=0}^{1-\epsilon}(r')} {E^{2-\epsilon}(r')}\,dr'$$

.

$B_\perp = |\vec{B}_\perp|$; $\vec{B}_\perp = \vec{B}_{wire}+\vec{B}_\mu$;

$\vec{B}_{wire}$: magnetic field parallel to the anode wire of the given tube;

$\vec{B}_\mu$: magnetic field magnetic field perpendicular to wire and parallel to the muon trajectory in the given tube.

Definition at line 47 of file BFieldCorFunc.h.

Member Typedef Documentation

ParVec

using MuonCalib::CalibFunc::ParVec = std::vector<double>

inherited

Definition at line 35 of file CalibFunc.h.

Constructor & Destructor Documentation

BFieldCorFunc() [1/2]

BFieldCorFunc::BFieldCorFunc ( const std::string &  quality, const CalibFunc::ParVec &  parameters, const IRtRelation *  rt  )

explicit

Constructor: quality = "high", slow but accurate initialization initialization of the correction function, quality = "medium", compromise between speed and accuracy of the initialization of the correction function (default), quality = "low", fast initialization of the correction function at the price of lower quality.

parameters[0] = high voltage [V], parameters[1] = $\epsilon$,

Definition at line 20 of file BFieldCorFunc.cxx.

                                                                                                               :
    IMdtBFieldCorFunc(parameters) {
    init(quality, parameters, rt);
}

BFieldCorFunc() [2/2]

BFieldCorFunc::BFieldCorFunc ( const CalibFunc::ParVec &  parameters, const IRtRelation *  rt  )

explicit

Constructor: parameters[0] = high voltage [V], parameters[1] = $\epsilon$,.

Definition at line 25 of file BFieldCorFunc.cxx.

                                                                                     : 
    IMdtBFieldCorFunc(parameters) {
    init("medium", parameters, rt);
}

Member Function Documentation

correction()

double BFieldCorFunc::correction ( double  t, double  B_wire, double  B_mu  ) const

virtual

get t(r, $\vec{B}$!=0)-t(r, $\vec{B}$=0); t = drift time t [ns] for B=0; B_wire = magnetic field parallel to the anode wire of the given tube, B_mu = magnetic field orthogonal to the wire and parallel to the muon trajectory in the given tube [B] = Tesla

Implements MuonCalib::IMdtBFieldCorFunc.

Definition at line 182 of file BFieldCorFunc.cxx.

                                                                           {
    // VARIABLES //
    double B_perp{std::hypot(B_wire, B_mu)};  // B orthogonal to the
                                              // electron drift path
    double B_factor{std::pow(B_perp, 2.0 - m_param[1])};
    double precision{0.1};                                             // precision of the correction in ns
    double t_max{t};                                                   // upper time search limit
    double t_min{t - 2 * correction_to_B(t, B_wire, B_mu, B_factor)};  // lower time search limit
    if (t_min < m_t_min) t_min = m_t_min;
    double time{t};                           // auxiliary time variable
    double integ{0.0};                        // integral
    double tmean{0.5 * (m_t_min + m_t_max)};  // mean time
    double tlength{m_t_max - m_t_min};        // length of drift-time interval
 
    // DRIFT TIME CHECK //
    if (t <= m_t_min) { return 0.0; }
    if (t > m_t_max) {
        t_max = m_t_max;
        time = m_t_max;
    }
 
    // SEARCH FOR THE CORRECTED DRIFT TIME //
    while (t_max - t_min > precision) {
        integ = 0.0;
        for (int k = 0; k < m_alpha.rows(); k++) {
            integ += m_alpha[k] * std::legendre(k, 2 * (0.5 * (t_min + t_max) - tmean) / tlength);
        }
        if (0.5 * (t_min + t_max) + B_factor * integ > time) {
            t_max = 0.5 * (t_min + t_max);
        } else {
            t_min = 0.5 * (t_min + t_max);
        }
    }
 
    return B_factor * integ;
}  // end BFieldCorFunc::correction

correction_to_B()

double BFieldCorFunc::correction_to_B	(	double	t,
		double	B_wire,
		double	B_mu,
		double	B_factor = -1.0
	)		const

Definition at line 225 of file BFieldCorFunc.cxx.

                                                                                                 {
    // VARIABLES //
    if (B_factor < 0) {
        const double B_perp{std::hypot(B_wire, B_mu)};
        // B orthogonal to the electron drift path
        B_factor = std::pow(B_perp, 2.0 - m_param[1]);
    }
    double time{t};
    double integ{0.0};                        // integral
    double tmean{0.5 * (m_t_min + m_t_max)};  // mean time
    double tlength{m_t_max - m_t_min};        // length of drift-time interval
 
    // DRIFT TIME CHECK //
    if (t <= m_t_min) { return 0.0; }
    if (t > m_t_max) { time = m_t_max; }
 
    // CALCULATE THE CORRECTION //
    for (int k = 0; k < m_alpha.rows(); k++) { 
        integ += m_alpha[k] * std::legendre(k, 2 * (time - tmean) / tlength); 
    }
 
    return B_factor * integ;
}  // end BFieldCorFunc::correction_to_B

epsilon()

double BFieldCorFunc::epsilon

(

)

const

< get the $\epsilon$ parameter of the B-field correction function

Definition at line 175 of file BFieldCorFunc.cxx.

{ return m_param[1]; }

init()

void BFieldCorFunc::init ( const std::string &  quality, const CalibFunc::ParVec &  params, const IRtRelation *  rt  )

private

Definition at line 30 of file BFieldCorFunc.cxx.

                                                                                                       {
    // PARAMETERS //
    m_quality = quality;
    m_param = params;
 
    // CONSISTENCY CHECK //
    if (m_param.size() != 2) {
        THROW_EXCEPTION("Wrong number of parameters!");
        return;
    }
 
    // VARIABLES //
    unsigned int nb_points(31);    // number of sample points for the integral
                                   // in the correction function
    double step{0.};                   // r step size
    double time{0.};                   // auxiliary time variable
    BaseFunctionFitter fitter(6);  // 6 fit parameters for the integral by
                                   // default ("medium quality")
    LegendrePolynomial legendre{};
 
    // QUALITY SETTING //
    if (m_quality == "high") {
        fitter.set_number_of_coefficients(8);
        nb_points = 31;
        m_step_size = 0.02;
    }
    if (m_quality == "medium") {
        fitter.set_number_of_coefficients(8);
        nb_points = 31;
        m_step_size = 0.06;
    }
    if (m_quality == "low") {
        fitter.set_number_of_coefficients(8);
        nb_points = 31;
        m_step_size = 0.12;
    }
    // sample points for the integral factor in the correction function
    std::vector<SamplePoint> sample_points(nb_points);
 
    // CALCULATE THE INTEGRAL PART OF THE CORRECTION FUNCTION //
    m_t_min = (rt)->tLower();
    m_t_max = (rt)->tUpper();
 
    // minimum and maximum radius //
    m_r_min = 0.025 * CLHEP::mm;    // minimum radius
    m_r_max = rt->radius(m_t_max);  // maximum radius
    if (m_r_max > 17.0 || m_r_max < m_r_min) {
        MsgStream log(Athena::getMessageSvc(), "BFieldCorFunc");
        log << MSG::INFO << "UNPHYSICAL MAXIMUM DRIFT RADIUS OF " << m_r_max << ", WILL BE SET TO 17.0!" << endmsg;
        m_r_max = 17.0;
    }
    step = ((m_r_max - m_r_min) / static_cast<double>(nb_points - 1));
 
    // set the sample points //
    double prev_r = 0;
    double prev_integral = 0;
    for (unsigned int k = 0; k < nb_points; k++) {
        time = t_from_r(m_r_min + k * step, rt);
        sample_points[k].set_x1(2 * (time - 0.5 * (m_t_min + m_t_max)) / (m_t_max - m_t_min));
        double new_r = rt->radius(time);
        double new_integral = 1.0e9 * integral(prev_r, new_r, rt) + prev_integral;
        sample_points[k].set_x2(new_integral);
        sample_points[k].set_error(1.0);
        prev_r = new_r;
        prev_integral = new_integral;
    }
 
    // perform the fit //
    fitter.fit_parameters(sample_points, 1, nb_points, legendre);
    m_alpha = fitter.coefficients();
 
}  // end BFieldCorFunc::init

integral()

double BFieldCorFunc::integral ( const double  r_min, const double  r_max, const IRtRelation *  rt  ) const

private

Definition at line 140 of file BFieldCorFunc.cxx.

                                                                                                  {
    // catch fp exceptions//
    if (m_r_min < 1e-10 || r_min < m_r_min) return 0.0;
 
    // VARIABLES //
    const double E0{m_param[0] / std::log(m_r_max / m_r_min)};  // E(r)=E0/r
    double radius{r_max};
    double rp{r_min};                                     // auxiliary radius variables
    double integ{0.0};                                    // current value of the integral
                                                          //    double step(0.010); // integration step size [mm]
    double step{m_step_size};                             // integration step size [mm]
    double time{0.};                                      // drift time
 
    // r IN [m_r_min, m_r_max]? //
    if (r_max < r_min) { return 0.0; }
    if (r_max > m_r_max) { radius = m_r_max; }
 
    // INTEGRATION //
    double delta = step;
    while (rp < radius) {
        time = t_from_r(rp, rt);
        if (rp + step > radius) delta = radius - rp;
        integ += 1.0e-3 * delta * std::pow(std::abs(rt->driftVelocity(time)) * 1.0e6, 1.0 - m_param[1]) /
                            std::pow(E0 / (rp * 1.0e-3), 2.0 - m_param[1]);
        rp += step;
    }
 
    return integ;
}  // end BFieldCorFunc::integral

name()

std::string BFieldCorFunc::name ( ) const

virtual

get the class name

get t(r, $\vec{B}$ !=0)-t(r, $\vec{B}$ =0); t = measured drift time t [ns]; B_wire = magnetic field parallel to the anode wire of the given tube, B_mu = magnetic field orthogonal to the wire and parallel to the muon trajectory in the given tube [B] = Tesla

Implements MuonCalib::CalibFunc.

Definition at line 180 of file BFieldCorFunc.cxx.

{ return std::string("BFieldCorFunc"); }

nPar()

unsigned int MuonCalib::CalibFunc::nPar ( ) const

inlineinherited

Definition at line 39 of file CalibFunc.h.

{ return m_parameters.size(); }

par()

double MuonCalib::CalibFunc::par ( unsigned int  index ) const

inlineinherited

Definition at line 41 of file CalibFunc.h.

                                                 {
                return index < nPar() ? m_parameters[index] : 0.;
            }

parameters()

const ParVec& MuonCalib::CalibFunc::parameters ( ) const

inlineinherited

Definition at line 40 of file CalibFunc.h.

{ return m_parameters; }

setRtRelationship()

void BFieldCorFunc::setRtRelationship

(

const IRtRelation &

rt

)

< set the $\epsilon$ parameter of the B-field correction function = eps

< set the r-t relationship used to calculate the B field correction to the measured drift time = rt

Definition at line 176 of file BFieldCorFunc.cxx.

                                                           {
    init(m_quality, m_param, &rt);
}

t_from_r()

double BFieldCorFunc::t_from_r ( const double  r, const IRtRelation *  rt  ) const

private

Definition at line 112 of file BFieldCorFunc.cxx.

                                                                          {
    // VARIABLES //
    constexpr double precision{0.010};  // spatial precision of the inversion
    double t_max{m_t_max};    // upper time search limit
    double t_min{m_t_min};    // lower time search limit
    double r_max{m_r_max};    // upper radius search limit
    double r_min{m_r_min};    // lower radius search limit
                              // SEARCH FOR THE CORRESPONDING DRIFT TIME //
    double t_guess, r_guess;
 
    do {
        t_guess = t_min + (t_max - t_min) / (r_max - r_min) * (r - r_min);
        r_guess = rt->radius(t_guess);
        if (r_guess > r) {
            r_max = r_guess;
            t_max = t_guess;
        } else {
            r_min = r_guess;
            t_min = t_guess;
        }
    } while (t_max - t_min > 0.1 && std::abs(r_guess - r) > precision);
    return t_guess;
}  // end BFieldCorFunc::t_from_r

typeName()

virtual std::string MuonCalib::IMdtBFieldCorFunc::typeName ( ) const

inlinevirtualinherited

Implements MuonCalib::CalibFunc.

Definition at line 18 of file IMdtBFieldCorFunc.h.

{ return "IMdtBFieldCorFunc"; }

Member Data Documentation

m_alpha

Amg::VectorX MuonCalib::BFieldCorFunc::m_alpha

private

Definition at line 109 of file BFieldCorFunc.h.

m_param

std::vector<double> MuonCalib::BFieldCorFunc::m_param {}

private

Definition at line 102 of file BFieldCorFunc.h.

m_parameters

ParVec MuonCalib::CalibFunc::m_parameters {}

privateinherited

Definition at line 48 of file CalibFunc.h.

m_quality

std::string MuonCalib::BFieldCorFunc::m_quality {}

private

Definition at line 105 of file BFieldCorFunc.h.

m_r_max

double MuonCalib::BFieldCorFunc::m_r_max {std::numeric_limits<double>::max()}

private

Definition at line 120 of file BFieldCorFunc.h.

m_r_min

double MuonCalib::BFieldCorFunc::m_r_min {-std::numeric_limits<double>::max()}

private

Definition at line 119 of file BFieldCorFunc.h.

m_step_size

double MuonCalib::BFieldCorFunc::m_step_size {0.}

private

Definition at line 106 of file BFieldCorFunc.h.

m_t_max

double MuonCalib::BFieldCorFunc::m_t_max {std::numeric_limits<double>::max()}

private

Definition at line 117 of file BFieldCorFunc.h.

m_t_min

double MuonCalib::BFieldCorFunc::m_t_min {-std::numeric_limits<double>::max()}

private

Definition at line 116 of file BFieldCorFunc.h.

---

The documentation for this class was generated from the following files:

• BFieldCorFunc.h
• BFieldCorFunc.cxx