#include <RtCalibrationCurved.h>

Inheritance diagram for MuonCalib::RtCalibrationCurved:

Public Types
using	MuonSegVec = std::vector<std::shared_ptr<MuonCalibSegment>>
using	MuonSegIt = MuonSegVec::iterator
using	MuonSegCit = MuonSegVec::const_iterator
using	MdtCalibOutputPtr = std::shared_ptr<IMdtCalibrationOutput>

Public Member Functions
	RtCalibrationCurved (const std::string &name)
	Default constructor: r-t accuracy is set to 0.5 mm.
	RtCalibrationCurved (const std::string &name, const double rt_accuracy, const unsigned int &func_type, const unsigned int &ord, const bool &fix_min, const bool &fix_max, const int &max_it, bool do_parabolic_extrapolation=false, bool do_smoothing=false, bool do_multilayer_rt_scale=false)
	Constructor.
	~RtCalibrationCurved ()
	Destructor.
double	reliability () const
	get the reliability of the final r-t relationship: 0: no convergence yet 1: convergence, r(t) is reliable 2: convergence, r(t) is unreliable
double	estimatedRtAccuracy () const
	get the estimated r-t quality (CLHEP::mm), the accuracy of the input r-t is computed at the end of the iteration; in order to get the accuracy of the final r-t, the algorithm has to be rerun with the final r-t as an input
int	numberOfSegments () const
	get the number of segments which were passed to the algorithm
int	numberOfSegmentsUsed () const
	get the number of segments which are used in the autocalibration
int	iteration () const
	get the number of the current iteration
bool	smoothing () const
	returns true, if the r-t relationship will be smoothened using the conventional autocalibration after convergence; returns false otherwise
void	setEstimateRtAccuracy (const double acc)
	set the estimated r-t accuracy =acc
void	fullMatrix (const bool &yes_or_no)
	yes_or_no=true: the full matrix relating the errors in the r-t relationship to the residuals is used yes_or_no=false: unit matrix is used (algorithm is equivalent to to the conventional/classical method
void	switch_on_control_histograms (const std::string &file_name)
	this methods requests control histograms from the algorithms; the algorithm will write them to ROOT file called "file_name"
void	switch_off_control_histograms ()
	the algorithm does not produce controll histograms (this is the default)
void	forceMonotony ()
	force r(t) to be monotonically increasing (this is default)
void	doNotForceMonotony ()
	do not force r(t) to be monotonically increasing
void	doParabolicExtrapolation ()
	requires that parabolic extrapolation will be used for small and large radii
void	noParabolicExtrapolation ()
	no parabolic extrapolation is done
void	doSmoothing ()
	requires that the r-t relationship will be smoothened using the conventional autocalibration after convergence
void	noSmoothing ()
	do not smoothen the r-t relationship after convergence
MdtCalibOutputPtr	analyseSegments (const MuonSegVec &seg) override
	perform the full autocalibration including iterations (required since MdtCalibInterfaces-00-01-06)
bool	handleSegment (MuonCalibSegment &seg)
	analyse the segment "seg" (this method was required before MdtCalibInterfaces-00-01-06)
void	setInput (const IMdtCalibrationOutput *rt_input) override
	set the r-t relationship, the internal autocalibration objects are reset
bool	analyse (const MuonSegVec &seg)
	perform the autocalibration with the segments acquired so far
bool	converged () const
	returns true, if the autocalibration has converged
virtual MdtCalibOutputPtr	getResults () const override
	returns the final r-t relationship
virtual std::string	name () const
	returns name (region) of instance

Private Member Functions
void	init (const double rt_accuracy, const unsigned int &func_type, const unsigned int &ord, const bool &fix_min, const bool &fix_max, const int &max_it, bool do_parabolic_extrapolation, bool do_smoothing, bool do_multilayer_rt_scale)
double	t_from_r (const double r)
void	display_segment (MuonCalibSegment segment, std::ofstream &outfile, const CurvedLine curved_segment)
std::shared_ptr< RtRelationLookUp >	performParabolicExtrapolation (const bool &min, const bool &max, const IRtRelation &in_rt)

Private Attributes
bool	m_control_histograms = false
bool	m_fix_min = false
bool	m_fix_max = false
int	m_max_it = 0
bool	m_force_monotony = false
bool	m_do_multilayer_rt_scale = false
int	m_nb_segments = 0
int	m_nb_segments_used = 0
int	m_iteration = 0
std::array< bool, 2 >	m_multilayer {}
int	m_status = 0
double	m_rt_accuracy = 0.0
double	m_rt_accuracy_previous
double	m_chi2_previous = 0.0
double	m_chi2 = 0.0
std::shared_ptr< const IRtRelation >	m_rt
double	m_t_length = 0.0
double	m_t_mean = 0.0
std::shared_ptr< IRtRelation >	m_rt_new
std::shared_ptr< RtCalibrationOutput >	m_output
std::unique_ptr< MultilayerRtDifference >	m_multilayer_rt_difference
double	m_r_max = 0.0
std::unique_ptr< CurvedPatRec >	m_tracker
CLHEP::HepSymMatrix	m_M_track
CLHEP::HepSymMatrix	m_M_track_inverse
bool	m_do_parabolic_extrapolation = false
bool	m_do_smoothing = false
unsigned int	m_order = 0U
std::vector< CLHEP::HepVector >	m_U
std::vector< CLHEP::HepVector >	m_U_weighted
CLHEP::HepSymMatrix	m_A
CLHEP::HepVector	m_alpha
CLHEP::HepVector	m_b
std::unique_ptr< BaseFunction >	m_base_function
std::unique_ptr< TFile >	m_tfile {}
std::unique_ptr< TH1F >	m_cut_evolution {}
std::unique_ptr< TH1F >	m_nb_segment_hits {}
std::unique_ptr< TH1F >	m_pull_initial {}
std::unique_ptr< TH1F >	m_pull_final {}
std::unique_ptr< TH2F >	m_residuals_initial {}
std::unique_ptr< TH2F >	m_residuals_initial_all {}
std::unique_ptr< TH2F >	m_residuals_final {}
std::unique_ptr< TH2F >	m_driftTime_initial {}
std::unique_ptr< TH2F >	m_driftTime_final {}
std::unique_ptr< TH2F >	m_adc_vs_residual_final {}
std::string	m_name

Detailed Description

Definition at line 45 of file RtCalibrationCurved.h.

Member Typedef Documentation

◆ MdtCalibOutputPtr

using MuonCalib::IMdtCalibration::MdtCalibOutputPtr = std::shared_ptr<IMdtCalibrationOutput>

inherited

Definition at line 30 of file IMdtCalibration.h.

◆ MuonSegCit

using MuonCalib::IMdtCalibration::MuonSegCit = MuonSegVec::const_iterator

inherited

Definition at line 29 of file IMdtCalibration.h.

◆ MuonSegIt

using MuonCalib::IMdtCalibration::MuonSegIt = MuonSegVec::iterator

inherited

Definition at line 28 of file IMdtCalibration.h.

◆ MuonSegVec

using MuonCalib::IMdtCalibration::MuonSegVec = std::vector<std::shared_ptr<MuonCalibSegment>>

inherited

Definition at line 27 of file IMdtCalibration.h.

Constructor & Destructor Documentation

◆ RtCalibrationCurved() [1/2]

RtCalibrationCurved::RtCalibrationCurved ( const std::string & name )

Default constructor: r-t accuracy is set to 0.5 mm.

The r-t accuracy is used internally to distinguish between good and bad segments. By default Legendre polynomials are used to parametrize the r-t correction. The order of the r-t correction polynomial is set to 15. Segments are not restricted to single multilayers. The full matrix relating the errors in r(t) to the residuals is used. No parabolic extrapolations are used. By default no smoothing is applied after convergence.

Definition at line 35 of file RtCalibrationCurved.cxx.

                                                              : IMdtCalibration(name), m_rt_accuracy_previous(0.0) {
    init(0.5 * CLHEP::mm, 1, 15, true, false, 15, false, false, false);
}

◆ RtCalibrationCurved() [2/2]

RtCalibrationCurved::RtCalibrationCurved	(	const std::string &	name,
		const double	rt_accuracy,
		const unsigned int &	func_type,
		const unsigned int &	ord,
		const bool &	fix_min,
		const bool &	fix_max,
		const int &	max_it,
		bool	do_parabolic_extrapolation = false,
		bool	do_smoothing = false,
		bool	do_multilayer_rt_scale = false )

Constructor.

Parameters

rt_accuracy	r-t accuracy in mm. The r-t accuracy is used internally to distinguish between good and bad segments.
func_type	Type of function to be used for the r-t correction; = 1: Legendre polynomial (default), = 2: Chebyshev polynomial, = 3: polygon equidistant in r.
ord	Order of the r-t correction function (15 is default).
fix_min	=true: r(t0) is fixed in the autocalibration procedure (this is default).
fix_max	=true: r(tmax) is fixed in the autocalibration procedure (false is default).
max_it	Maximum number of iterations (20 by default).
do_parabolic_extrapolation	Use parabolic extrapolations for small and larged drift radii.
do_smoothing	Smoothing of the r-t relations after convergence.

Definition at line 39 of file RtCalibrationCurved.cxx.

                                                                                                                          :
    IMdtCalibration(name), m_rt_accuracy_previous(0.0) {
    init(rt_accuracy, func_type, ord, fix_min, fix_max, max_it, do_parabolic_extrapolation, do_smoothing, do_multilayer_rt_scale);
}

◆ ~RtCalibrationCurved()

RtCalibrationCurved::~RtCalibrationCurved ( )

Destructor.

Definition at line 46 of file RtCalibrationCurved.cxx.

                                          {
    if (m_tfile) { m_tfile->Write(); }
}

Member Function Documentation

◆ analyse()

bool RtCalibrationCurved::analyse ( const MuonSegVec & seg )

perform the autocalibration with the segments acquired so far

Definition at line 563 of file RtCalibrationCurved.cxx.

                                                       {
    // VARIABLES //
    int ifail;                   // flag indicating a failure of the matrix inversion
    unsigned int nb_points(30);  // number of points used to set the new r-t relationship
    double step;                 // r step size
    auto rt_Chebyshev = dynamic_cast<const RtChebyshev*>(m_rt.get());
    auto rt_LookUp = dynamic_cast<const RtRelationLookUp*>(m_rt.get());
    double r_corr;                               // radial correction
    std::vector<double> rt_param(m_rt->nPar());  // parameters for the new r-t
    double x;                                    // reduced time
    // SOLVE THE AUTOCALIBRATION EQUATION //
 
    m_alpha = m_A.inverse(ifail) * m_b;
    if (ifail != 0) {
        MsgStream log(Athena::getMessageSvc(), "RtCalibrationCurved");
        log << MSG::WARNING << "analyse() - Could not solve the autocalibration equation!" << endmsg;
        return false;
    }
    // CALCULATE THE NEW r-t RELATIONSHIP //
 
    // input r-t is of type RtChebyshev //
    if (rt_Chebyshev) {
        // set the number of points //
        if (rt_Chebyshev->nDoF() > 30) { nb_points = rt_Chebyshev->nDoF(); }
 
        // r step size //
        step = m_r_max / static_cast<double>(nb_points);
 
        // sample points and Chebyshev fitter //
        std::vector<SamplePoint> x_r(nb_points + 1);
        BaseFunctionFitter fitter(rt_Chebyshev->nDoF());
        ChebyshevPolynomial chebyshev;
 
        // calculate the sample points //
        for (unsigned int k = 0; k < nb_points + 1; k++) {
            x_r[k].set_x1(t_from_r(k * step));
            x_r[k].set_x2(rt_Chebyshev->radius(x_r[k].x1()));
            x_r[k].set_x1(rt_Chebyshev->getReducedTime(x_r[k].x1()));
            x_r[k].set_error(1.0);
 
            r_corr = 0.0;
            for (unsigned int l = 0; l < m_order; l++) {
                //              r_corr = r_corr+m_alpha[l]*
                //                  m_base_function->value(l, x_r[k].x1());
                r_corr = r_corr + m_alpha[l] * m_base_function->value(l, (x_r[k].x2() - 0.5 * m_r_max) / (0.5 * m_r_max));
            }
 
            // do not change the r-t and the endpoints //
            if (((k == 0 || x_r[k].x2() < 0.5) && m_fix_min) || ((k == nb_points || x_r[k].x2() > m_r_max) && m_fix_max)) {
                r_corr = 0.0;
                x_r[k].set_error(0.01);
            }
            x_r[k].set_x2(x_r[k].x2() - r_corr);
        }
 
        // force monotony //
        if (m_force_monotony) {
            for (unsigned int k = 0; k < nb_points; k++) {
                if (x_r[k].x2() > x_r[k + 1].x2()) { x_r[k + 1].set_x2(x_r[k].x2()); }
            }
        }
 
        // create the new r-t relationship //
        fitter.fit_parameters(x_r, 1, nb_points + 1, chebyshev);
        rt_param[0] = rt_Chebyshev->tLower();
        rt_param[1] = rt_Chebyshev->tUpper();
        for (unsigned int k = 0; k < rt_Chebyshev->nDoF(); k++) { rt_param[k + 2] = fitter.coefficients()[k]; }
 
        m_rt_new = std::make_shared<RtChebyshev>(rt_param);
        m_rt_new->SetTmaxDiff(m_rt->GetTmaxDiff());
    }
 
    // input-rt is of type RtRelationLookUp //
    if (rt_LookUp) {
        rt_param = rt_LookUp->parameters();
        unsigned int min_k(2), max_k(rt_param.size());
        if (m_fix_min) { min_k = 3; }
        if (m_fix_max) { max_k = rt_param.size() - 1; }
        for (unsigned int k = min_k; k < max_k; k++) {
            x = (rt_param[k] - 0.5 * m_r_max) / (0.5 * m_r_max);
            r_corr = m_alpha[0];
            for (unsigned int l = 1; l < m_order; l++) { r_corr = r_corr + m_alpha[l] * m_base_function->value(l, x); }
            rt_param[k] = rt_param[k] - r_corr;
            if (m_force_monotony && k > 2) {
                if (rt_param[k] < rt_param[k - 1]) { rt_param[k] = rt_param[k - 1]; }
            }
        }
 
        m_rt_new = std::make_shared<RtRelationLookUp>(rt_param);
        m_rt_new->SetTmaxDiff(m_rt->GetTmaxDiff());
    }
    // DETERMINE THE r-t QUALITY AND CHECK FOR CONVERGENCE //
 
    // estimate r-t accuracy //
    m_rt_accuracy = 0.0;
    //  double m_rt_accuracy_diff = 0.0;
    double r_corr_max = 0.0;
 
    for (unsigned int k = 0; k < 100; k++) {
        r_corr = m_alpha[0];
        for (unsigned int l = 1; l < m_order; l++) { r_corr = r_corr + m_alpha[l] * m_base_function->value(l, -1.0 + 0.01 * k); }
        if (std::abs(r_corr) > r_corr_max) { r_corr_max = std::abs(r_corr); }
        m_rt_accuracy = m_rt_accuracy + r_corr * r_corr;
    }
    m_rt_accuracy = std::sqrt(0.01 * m_rt_accuracy);
    //  m_rt_accuracy_diff = m_rt_accuracy_previous - m_rt_accuracy;
    m_rt_accuracy_previous = m_rt_accuracy;
 
    // convergence? //
    m_chi2 = m_chi2 / static_cast<double>(m_nb_segments_used);
    if ((m_chi2 <= m_chi2_previous || std::abs(m_chi2 - m_chi2_previous) > 0.01) ||
        (std::abs(m_rt_accuracy) > 0.001 && m_iteration < m_max_it)) {
        m_status = 0;  // no convergence yet
    } else {
        m_status = 1;
    }
    m_chi2_previous = m_chi2;
 
    // reliabilty of convergence //
    if (m_status != 0) {
        if (m_nb_segments_used < 0.5 * m_nb_segments) { m_status = 2; }
    }
    // Set new multilayer rt-difference //
    if (m_iteration > 0 && m_do_multilayer_rt_scale) {
        m_multilayer_rt_difference->DoFit(m_rt_new.get(), seg);
    } else {
        m_multilayer_rt_difference->DoFit(nullptr, {});
    }
    // STORE THE RESULTS IN THE RtCalibrationOutput //
 
    m_iteration = m_iteration + 1;
 
    m_output = std::make_shared<RtCalibrationOutput>(
        m_rt_new, std::make_shared<RtFullInfo>("RtCalibrationCurved", m_iteration, m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
 
    return true;
}

◆ analyseSegments()

RtCalibrationCurved::MdtCalibOutputPtr RtCalibrationCurved::analyseSegments ( const MuonSegVec & seg )

overridevirtual

perform the full autocalibration including iterations (required since MdtCalibInterfaces-00-01-06)

Implements MuonCalib::IMdtCalibration.

Definition at line 101 of file RtCalibrationCurved.cxx.

                                                                                               {
    std::shared_ptr<const IRtRelation> tmp_rt;
    // Initial RESIDUALS //
    for (const auto & k : seg) {
        for (unsigned int l = 0; l < k->hitsOnTrack(); l++) {
            double rr = (k->mdtHOT())[l]->driftRadius();
            double dd = (k->mdtHOT())[l]->signedDistanceToTrack();
            if (m_residuals_final) m_residuals_initial_all->Fill(std::abs(dd), std::abs(rr) - std::abs(dd));
            m_driftTime_initial->Fill(rr, dd);
        }
    }
    // AUTOCALIBRATION LOOP //
    while (!converged()) {
        for (const auto & k : seg) { handleSegment(*k); }
        if (!analyse(seg)) {
            MsgStream log(Athena::getMessageSvc(), "RtCalibrationCurved");
            log << MSG::WARNING << "analyseSegments() - analyse failed, segments:" << endmsg;
            for (unsigned int i = 0; i < seg.size(); i++) {
                log << MSG::WARNING << i << " " << seg[i]->direction() << " " << seg[i]->position() << endmsg;
            }
            return nullptr;
        }
 
        const RtCalibrationOutput *rtOut = m_output.get();
 
        if (!converged()) { setInput(rtOut); }
        tmp_rt = rtOut->rt();
    }
 
    // parabolic extrapolations for small radii //
    if (m_do_parabolic_extrapolation) {
        std::shared_ptr<RtRelationLookUp> tmprt = performParabolicExtrapolation(true, true, *tmp_rt);
        m_output = std::make_shared<RtCalibrationOutput>(
            tmprt, std::make_shared<RtFullInfo>("RtCalibrationCurved", m_iteration, m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
 
        tmp_rt = tmprt;
    }
    // SMOOTHING AFTER CONVERGENCE IF REQUESTED //
    if (!m_do_smoothing) {
        // final residuals //
        double r{0}, d{0}, adc{0};
        for (const auto & k : seg) {
            for (unsigned int l = 0; l < k->hitsOnTrack(); l++) {
                adc = (k->mdtHOT())[l]->adcCount();
                r = (k->mdtHOT())[l]->driftRadius();
                d = (k->mdtHOT())[l]->signedDistanceToTrack();
                if (m_residuals_final) m_residuals_final->Fill(std::abs(d), std::abs(r) - std::abs(d));
                m_driftTime_final->Fill(r, d);
                m_adc_vs_residual_final->Fill(adc, std::abs(r) - std::abs(d));
            }
        }
        return getResults();
    }
 
    // maximum number of iterations //
    int max_smoothing_iterations(static_cast<int>(m_max_it));
    if (max_smoothing_iterations == 0) { max_smoothing_iterations = 1; }
 
    // convergence RMS //
    double convergence_RMS(0.002);
 
    // variables //
    int it(0);        // iteration
    double RMS(1.0);  // RMS change of r(t)
 
    // smoothing //
    //---------------------------------------------------------------------------//
    //---------------------------------------------------------------------------//
    while (it < max_smoothing_iterations && RMS > convergence_RMS) {
        //---------------------------------------------------------------------------//
        AdaptiveResidualSmoothing smoothing;
 
        // counter //
        unsigned int counter{0};
        // overwrite drift radii and calculate the average resolution //
        for (const auto & k : seg) {
            if (k->mdtHitsOnTrack() < 4) { continue; }
            double avres(0.0);
            for (unsigned int h = 0; h < k->mdtHitsOnTrack(); h++) {
                k->mdtHOT()[h]->setDriftRadius(tmp_rt->radius(k->mdtHOT()[h]->driftTime()),
                                                    k->mdtHOT()[h]->sigmaDriftRadius());
                if (k->mdtHOT()[h]->sigmaDriftRadius() < 0.5 * m_r_max) {
                    avres = avres + k->mdtHOT()[h]->sigma2DriftRadius();
                } else {
                    avres = avres + 0.01;
                }
            }
            avres = avres / static_cast<double>(k->mdtHitsOnTrack());
            avres = std::sqrt(avres);
            if (smoothing.addResidualsFromSegment(*k, true, 5.0 * avres)) { counter++; }
        }
 
        // break, do no smoothing if there are not enough segments //
        if (counter < 1000) {
            MsgStream log(Athena::getMessageSvc(), "RtCalibrationCurved");
            log << MSG::WARNING << "analyseSegments() - too small number of reconstructed segments!" << endmsg;
            // final residuals //
            for (const auto & k : seg) {
                for (unsigned int l = 0; l < k->hitsOnTrack(); l++) {
                    double r = (k->mdtHOT())[l]->driftRadius();
                    double d = (k->mdtHOT())[l]->signedDistanceToTrack();
                    if (m_residuals_final) m_residuals_final->Fill(d, std::abs(r) - std::abs(d), 1.0);
                }
            }
            return getResults();
        }
 
        // smoothing //
        RtRelationLookUp smooth_rt(smoothing.performSmoothing(*tmp_rt, m_fix_min, m_fix_max));
 
        // calculate RMS change //
        RMS = 0.0;
        double bin_width(0.01 * (smooth_rt.tUpper() - smooth_rt.tLower()));
        for (double t = smooth_rt.tLower(); t <= smooth_rt.tUpper(); t = t + bin_width) {
            RMS = RMS + std::pow(smooth_rt.radius(t) - tmp_rt->radius(t), 2);
        }
        RMS = std::sqrt(0.01 * RMS);
 
        // increase the iterations counter //
        it++;
 
        // delete tmp_rt and update it //
        tmp_rt = std::make_shared<RtRelationLookUp>(smooth_rt);
 
        //---------------------------------------------------------------------------//
    }
    //---------------------------------------------------------------------------//
    //---------------------------------------------------------------------------//
    m_output = std::make_shared<RtCalibrationOutput>(
        tmp_rt, std::make_shared<RtFullInfo>("RtCalibrationCurved", m_iteration, m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
    // FINAL RESIDUALS //
    for (const auto & k : seg) {
        for (unsigned int l = 0; l < k->hitsOnTrack(); l++) {
            double adc = (k->mdtHOT())[l]->adcCount();
            double r = (k->mdtHOT())[l]->driftRadius();
            double d = (k->mdtHOT())[l]->signedDistanceToTrack();
            if (m_residuals_final) m_residuals_final->Fill(d, std::abs(r) - std::abs(d));
            m_driftTime_final->Fill(r, d);
            m_adc_vs_residual_final->Fill(adc, std::abs(r) - std::abs(d));
        }
    }
    // RETURN THE RESULT AFTER CONVERGENCE //
    return getResults();
}

◆ converged()

bool RtCalibrationCurved::converged ( ) const

returns true, if the autocalibration has converged

Definition at line 717 of file RtCalibrationCurved.cxx.

717{ return (m_status > 0); }

◆ display_segment()

void RtCalibrationCurved::display_segment	(	MuonCalibSegment *	segment,
		std::ofstream &	outfile,
		const CurvedLine *	curved_segment )

private

Definition at line 824 of file RtCalibrationCurved.cxx.

                                                                                                                           {
    // VARIABLES //
    double y_min, y_max, z_min, z_max;  // minimum and maximum y and z coordinates
    Amg::Vector3D null(0.0, 0.0, 0.0);  // auxiliary null vector
    // DISPLAY THE SEGMENT //
    // minimum and maximum coordinates //
    y_min = (segment->mdtHOT()[0])->localPosition().y();
    y_max = y_min;
    z_min = (segment->mdtHOT()[0])->localPosition().z();
    z_max = z_min;
    for (unsigned int k = 1; k < segment->mdtHitsOnTrack(); k++) {
        if ((segment->mdtHOT()[k])->localPosition().y() < y_min) { y_min = (segment->mdtHOT()[k])->localPosition().y(); }
        if ((segment->mdtHOT()[k])->localPosition().y() > y_max) { y_max = (segment->mdtHOT()[k])->localPosition().y(); }
        if ((segment->mdtHOT()[k])->localPosition().z() < z_min) { z_min = (segment->mdtHOT()[k])->localPosition().z(); }
        if ((segment->mdtHOT()[k])->localPosition().z() > z_max) { z_max = (segment->mdtHOT()[k])->localPosition().z(); }
    }
    for (unsigned int k = 0; k < segment->mdtCloseHits(); k++) {
        if ((segment->mdtClose()[k])->localPosition().y() < y_min) { y_min = (segment->mdtClose()[k])->localPosition().y(); }
        if ((segment->mdtClose()[k])->localPosition().y() > y_max) { y_max = (segment->mdtClose()[k])->localPosition().y(); }
        if ((segment->mdtClose()[k])->localPosition().z() < z_min) { z_min = (segment->mdtClose()[k])->localPosition().z(); }
        if ((segment->mdtClose()[k])->localPosition().z() > z_max) { z_max = (segment->mdtClose()[k])->localPosition().z(); }
    }
 
    // write out the coordinate system //
    if (y_max - y_min > z_max - z_min) {
        outfile << "nullptr " << y_min - 30.0 << " " << y_max + 30.0 << " " << 0.5 * (z_min + z_max) - 0.5 * (y_max - y_min) - 30.0 << " "
                << 0.5 * (z_min + z_max) + 0.5 * (y_max - y_min) + 30.0 << "\n";
    } else {
        outfile << "nullptr " << 0.5 * (y_min + y_max) - 0.5 * (z_max - z_min) - 30.0 << " "
                << 0.5 * (y_min + y_max) + 0.5 * (z_max - z_min) + 30.0 << " " << z_min - 30.0 << " " << z_max + 30.0 << "\n";
    }
 
    // write out the hits on track //
    for (unsigned int k = 0; k < segment->mdtHitsOnTrack(); k++) {
        // draw a circle for the tube //
        outfile << "SET PLCI 1\n"
                << "ARC " << (segment->mdtHOT()[k])->localPosition().y() << " " << (segment->mdtHOT()[k])->localPosition().z() << " 15.0\n";
 
        // draw a drift circle //
        outfile << "SET PLCI 3\n"
                << "ARC " << (segment->mdtHOT()[k])->localPosition().y() << " " << (segment->mdtHOT()[k])->localPosition().z() << " "
                << (segment->mdtHOT()[k])->driftRadius() << "\n";
    }
 
    // write out the close hits //
    for (unsigned int k = 0; k < segment->mdtCloseHits(); k++) {
        // draw a circle for the tube //
        outfile << "SET PLCI 1\n"
                << "ARC " << (segment->mdtClose()[k])->localPosition().y() << " " << (segment->mdtClose()[k])->localPosition().z()
                << " 15.0\n";
 
        // draw a drift circle //
        outfile << "SET PLCI 2\n"
                << "ARC " << (segment->mdtClose()[k])->localPosition().y() << " " << (segment->mdtClose()[k])->localPosition().z() << " "
                << (segment->mdtClose()[k])->driftRadius() << "\n";
    }
 
    // write out the reconstructed track //
    // a straight line is drawn by default //
    if (curved_segment == nullptr) {
        MTStraightLine aux_track(segment->position(), segment->direction(), null, null);
        outfile << "SET PLCI 4\n"
                << "LINE " << aux_track.a_x2() * (z_min - 30.0) + aux_track.b_x2() << " " << z_min - 30.0 << " "
                << aux_track.a_x2() * (z_max + 30.0) + aux_track.b_x2() << " " << z_max + 30.0 << "\n";
    }
 
    // a curved segment is drawn on demand //
    if (curved_segment != nullptr) {
        double step_size((60.0 + z_max - z_min) / 50.0);
        for (double aux_z = z_min; aux_z <= z_max; aux_z = aux_z + step_size) {
            outfile << "SET PLCI 4\n"
                    << "LINE " << curved_segment->getPointOnLine(aux_z).y() << " " << aux_z << " "
                    << curved_segment->getPointOnLine(aux_z + step_size).y() << " " << aux_z + step_size << "\n";
        }
    }
 
    // add a wait statement //
    outfile << "WAIT\n";
 
    return;
}

◆ doNotForceMonotony()

void RtCalibrationCurved::doNotForceMonotony ( )

do not force r(t) to be monotonically increasing

Definition at line 97 of file RtCalibrationCurved.cxx.

97{ m_force_monotony = false; }

◆ doParabolicExtrapolation()

void RtCalibrationCurved::doParabolicExtrapolation ( )

requires that parabolic extrapolation will be used for small and large radii

Definition at line 98 of file RtCalibrationCurved.cxx.

98{ m_do_parabolic_extrapolation = true; }

◆ doSmoothing()

void RtCalibrationCurved::doSmoothing ( )

requires that the r-t relationship will be smoothened using the conventional autocalibration after convergence

Definition at line 99 of file RtCalibrationCurved.cxx.

99{ m_do_smoothing = true; }

◆ estimatedRtAccuracy()

double RtCalibrationCurved::estimatedRtAccuracy ( ) const

get the estimated r-t quality (CLHEP::mm), the accuracy of the input r-t is computed at the end of the iteration; in order to get the accuracy of the final r-t, the algorithm has to be rerun with the final r-t as an input

Definition at line 51 of file RtCalibrationCurved.cxx.

51{ return m_rt_accuracy; }

◆ forceMonotony()

void RtCalibrationCurved::forceMonotony ( )

force r(t) to be monotonically increasing (this is default)

Definition at line 96 of file RtCalibrationCurved.cxx.

96{ m_force_monotony = true; }

◆ fullMatrix()

void MuonCalib::RtCalibrationCurved::fullMatrix ( const bool & yes_or_no )

yes_or_no=true: the full matrix relating the errors in the r-t relationship to the residuals is used yes_or_no=false: unit matrix is used (algorithm is equivalent to to the conventional/classical method

◆ getResults()

RtCalibrationCurved::MdtCalibOutputPtr RtCalibrationCurved::getResults ( ) const

overridevirtual

returns the final r-t relationship

Implements MuonCalib::IMdtCalibration.

Definition at line 718 of file RtCalibrationCurved.cxx.

718{ return m_output; }

◆ handleSegment()

bool RtCalibrationCurved::handleSegment ( MuonCalibSegment & seg )

analyse the segment "seg" (this method was required before MdtCalibInterfaces-00-01-06)

Definition at line 267 of file RtCalibrationCurved.cxx.

                                                             {
    // RETURN, IF THE SEGMENT HAD LESS THAN 4 HITS //
    if (m_control_histograms) { m_cut_evolution->Fill(0.0, 1.0); }
 
    m_nb_segments = m_nb_segments + 1;
    if (seg.mdtHitsOnTrack() < 4) { return true; }
 
    if (m_control_histograms) { m_cut_evolution->Fill(1.0, 1.0); }
 
    if (std::isnan(seg.direction().x()) || std::isnan(seg.direction().y()) || std::isnan(seg.direction().z()) ||
        std::isnan(seg.position().x()) || std::isnan(seg.position().y()) || std::isnan(seg.position().z())) {
        return true;
    }
    if (std::abs(seg.direction().y()) > 100) { return true; }
    // VARIABLES //
 
    // segment reconstruction and segment selection //
    double aux_res;            // auxiliary resolution
    double av_res(0.0);        // average spatial resolution of the tube hits
    double chi2_scale_factor;  // chi^2 scale factor used to take into
    // account bad r-t accuracy in the segment
    // selection
    IMdtSegmentFitter::HitSelection hit_selection(seg.mdtHitsOnTrack());
    // hit selection vectors for refits
    unsigned int nb_hits_in_ml;   // number of hits in the multilayers
    double x;                     // reduced time = (r(t)-0.5*m_r_max)/(0.5*m_r_max)
    std::vector<double> d_track;  // signed distances of the track from the anode
    // wires of the tubes
    std::vector<MTStraightLine> w;      // anode wires
    Amg::Vector3D null(0.0, 0.0, 0.0);  // auxiliary 0 vector
    Amg::Vector3D xhat(1.0, 0.0, 0.0);  // auxiliary unit vector
 
    // objects needed to calculate the autocalibration matrix and vector //
    std::vector<CLHEP::HepVector> F;  // auxiliary vectors for the calculation of the
    // track cooeffients matrix
    CLHEP::HepVector residual_value;     // residuals values
    CLHEP::HepVector weighted_residual;  // residual values weighted by the inverse
    // variance of the radius measurements
    CLHEP::HepMatrix D;  // Jacobian of the residuals (dresiduals/dr)
                         //     static ofstream display("display.kumac");
    // PREPARATION FOR THE SEGMENT REFIT //
 
    // overwrite the drift radii according to the input r-t relationship, //
    // calculate the average spatial resolution to define a road width,   //
    // select the hits in the chamber                                     //
    nb_hits_in_ml = 0;
    for (unsigned int k = 0; k < seg.mdtHitsOnTrack(); k++) {
        // make the resolution at small radii large enough //
        aux_res = seg.mdtHOT()[k]->sigmaDriftRadius();
        if (m_rt->radius(seg.mdtHOT()[k]->driftTime()) < 0.75) {
            if (aux_res < 0.5 - m_rt->radius(seg.mdtHOT()[k]->driftTime())) { aux_res = 0.5 - m_rt->radius(seg.mdtHOT()[k]->driftTime()); }
        }
 
        // overwrite radius //
        seg.mdtHOT()[k]->setDriftRadius(m_rt->radius(seg.mdtHOT()[k]->driftTime()), aux_res);
 
        // hit selection //
        hit_selection[k] = 0;
 
        // reject hits with bad radii //
        if (seg.mdtHOT()[k]->driftRadius() < 0.0 || seg.mdtHOT()[k]->driftRadius() > m_r_max ||
            seg.mdtHOT()[k]->sigmaDriftRadius() > 10.0 * m_r_max) {
            hit_selection[k] = 1;
        }
 
        // average resolution in the segment //
        if (hit_selection[k] == 0) {
            if (seg.mdtHOT()[k]->sigmaDriftRadius() < 0.5 * m_r_max) {
                av_res = av_res + std::pow(seg.mdtHOT()[k]->sigmaDriftRadius(), 2);
            } else {
                av_res = av_res + 0.01;
            }
        }
 
        // counting of selected segments //
        nb_hits_in_ml = nb_hits_in_ml + (1 - hit_selection[k]);
    }
 
    // average resolution and chi^2 scale factor //
    av_res = std::sqrt(av_res / static_cast<double>(seg.mdtHitsOnTrack()));
    chi2_scale_factor = std::hypot(av_res, m_rt_accuracy) / av_res;
    // FILL THE AUTOCALIBRATION MATRICES //
 
    // set the road width for the track reconstruction //
    m_tracker->setRoadWidth(7.0 * std::hypot(av_res, m_rt_accuracy));
 
    // check whether there are enough hits in the chambers //
    if (nb_hits_in_ml < 4) { return true; }
 
    // refit the segments //
    CurvedLine track;
 
    if (!m_tracker->fit(seg, hit_selection, track)) { return true; }
 
    if (std::isnan(track.chi2())) { return true; }
 
    // check the quality of the fit //
    if (track.chi2PerDegreesOfFreedom() > 5 * chi2_scale_factor) { return true; }
 
    // check whether there are at least three track hits //
    if (m_control_histograms) { m_nb_segment_hits->Fill(track.numberOfTrackHits()); }
    if (track.numberOfTrackHits() < 4) { return true; }
 
    //  display_segment(&seg, display&(m_tracker->curvedTrack()));
 
    // reject tracks with silly parameters
    if (std::abs(track.getTangent(seg.mdtHOT()[0]->localPosition().z()).a_x2()) > 8.0e8) { return true; }
 
    // bookkeeping for convergence decision and reliability estimation //
    m_chi2 += track.chi2PerDegreesOfFreedom();
    m_nb_segments_used = m_nb_segments_used + 1;
 
    // fill the autocalibration objects //
 
    // track coeffient matrix and its inverse //
    F = std::vector<CLHEP::HepVector>(track.numberOfTrackHits());
    for (unsigned int h = 0; h < track.numberOfTrackHits(); h++) {
        const MdtCalibHitBase &hb = *(track.trackHits()[h]);
        m_multilayer_rt_difference->Fill(hb, *m_rt);
 
        F[h] = CLHEP::HepVector(m_M_track.num_row());
        for (int p = 0; p < F[h].num_row(); p++) {
            double x = std::sqrt(1.0 + std::pow(track.getTangent((track.trackHits()[h]->localPosition()).z()).a_x2(), 2));
            if (x) {
                (F[h])[p] = std::legendre(p, (track.trackHits()[h]->localPosition()).z()) / x;
            } else {
                (F[h])[p] = 0.;
            }
        }
    }
    for (int p = 0; p < m_M_track.num_row(); p++) {
        for (int q = p; q < m_M_track.num_row(); q++) {
            m_M_track[p][q] = 0.0;
            for (unsigned int h = 0; h < track.numberOfTrackHits(); h++) {
                m_M_track[p][q] = m_M_track[p][q] + (F[h])[p] * (F[h])[q] / (track.trackHits()[h])->sigma2DriftRadius();
            }
            m_M_track[q][p] = m_M_track[p][q];
        }
    }
    int ifail;
    m_M_track_inverse = m_M_track.inverse(ifail);
    if (ifail != 0) {
        MsgStream log(Athena::getMessageSvc(), "RtCalibrationCurved");
        log << MSG::WARNING << "handleSegment() - Could not invert track matrix!" << endmsg;
        return true;
    }
 
    // Jacobian matrix for the residuals //
    // track distances, residuals, corrections //
    d_track = std::vector<double>(track.numberOfTrackHits());
    w = std::vector<MTStraightLine>(track.numberOfTrackHits());
    residual_value = CLHEP::HepVector(track.numberOfTrackHits());
    weighted_residual = CLHEP::HepVector(track.numberOfTrackHits());
    for (unsigned int l = 0; l < m_order; l++) {
        m_U[l] = CLHEP::HepVector(track.numberOfTrackHits());
        m_U_weighted[l] = CLHEP::HepVector(track.numberOfTrackHits());
    }
    for (unsigned int h = 0; h < track.numberOfTrackHits(); h++) {
        w[h] = MTStraightLine(Amg::Vector3D(0.0, (track.trackHits()[h]->localPosition()).y(), (track.trackHits()[h]->localPosition()).z()),
                              xhat, null, null);
        d_track[h] = track.getTangent((track.trackHits()[h]->localPosition()).z()).signDistFrom(w[h]);
        residual_value[h] = (track.trackHits()[h])->driftRadius() - std::abs(d_track[h]);
        if (m_control_histograms) {
            if (m_iteration == 0) { m_residuals_initial->Fill(std::abs(d_track[h]), residual_value[h], 1.0); }
        }
        for (unsigned int l = 0; l < m_order; l++) {
            x = (track.trackHits()[h]->driftRadius() - 0.5 * m_r_max) / (0.5 * m_r_max);
            (m_U[l])[h] = m_base_function->value(l, x);
        }
    }
 
    // Jacobian //
    D = CLHEP::HepMatrix(track.numberOfTrackHits(), track.numberOfTrackHits());
    for (unsigned int h = 0; h < track.numberOfTrackHits(); h++) {
        for (unsigned int hp = 0; hp < track.numberOfTrackHits(); hp++) {
            D[h][hp] = (h == hp) - (2 * (d_track[h] >= 0) - 1) * (2 * (d_track[hp] >= 0) - 1) * dot(F[h], m_M_track_inverse * F[hp]) /
                                       (track.trackHits()[hp])->sigma2DriftRadius();
        }
    }
 
    // autocalibration objects //
    // errors of the residuals //
    std::vector<double> sigma_residual(track.numberOfTrackHits(), 0.0);
    for (unsigned int h = 0; h < track.numberOfTrackHits(); h++) {
        for (unsigned int hp = 0; hp < track.numberOfTrackHits(); hp++) {
            sigma_residual[h] = sigma_residual[h] + std::pow(D[h][hp] * (track.trackHits()[hp])->sigmaDriftRadius(), 2);
        }
        sigma_residual[h] = std::sqrt(sigma_residual[h]);
        if (sigma_residual[h] < av_res / track.numberOfTrackHits()) { sigma_residual[h] = av_res / std::sqrt(track.numberOfTrackHits()); }
    }
 
    // weight residuals and correction functions //
    for (unsigned int h = 0; h < track.numberOfTrackHits(); h++) {
        weighted_residual[h] = residual_value[h] / sigma_residual[h];
        if (m_control_histograms) {
            if (m_iteration == 0) { m_pull_initial->Fill(weighted_residual[h], 1.0); }
        }
        for (unsigned int l = 0; l < m_order; l++) { (m_U_weighted[l])[h] = (m_U[l])[h] / sigma_residual[h]; }
    }
 
    // autocalibration matrix and autocalibration vector //
    CLHEP::HepSymMatrix A_tmp(m_A);
 
    // autocalibration objects //
    for (unsigned int p = 0; p < m_order; p++) {
        for (unsigned int pp = p; pp < m_order; pp++) {
            A_tmp[p][pp] = m_A[p][pp] + dot(D * m_U_weighted[p], D * m_U_weighted[pp]);
            if (std::isnan(A_tmp[p][pp])) { return true; }
        }
    }
 
    CLHEP::HepVector b_tmp(m_b);
    for (unsigned int p = 0; p < m_order; p++) {
        b_tmp[p] = m_b[p] + dot(D * m_U_weighted[p], weighted_residual);
        if (std::isnan(b_tmp[p])) { return true; }
    }
 
    m_A = A_tmp;
    m_b = b_tmp;
 
    return true;
}

◆ init()

void RtCalibrationCurved::init	(	const double	rt_accuracy,
		const unsigned int &	func_type,
		const unsigned int &	ord,
		const bool &	fix_min,
		const bool &	fix_max,
		const int &	max_it,
		bool	do_parabolic_extrapolation,
		bool	do_smoothing,
		bool	do_multilayer_rt_scale )

private

Definition at line 725 of file RtCalibrationCurved.cxx.

                                                            {
    // RESET PRIVATE VARIABLES //
    m_rt = nullptr;
    m_r_max = 15.0 * CLHEP::mm;
    m_control_histograms = false;
    m_nb_segments = 0;
    m_nb_segments_used = 0;
    m_iteration = 0;
    m_multilayer[0] = false;
    m_multilayer[1] = false;
    m_status = 0;
    m_rt_accuracy = std::abs(rt_accuracy);
    m_chi2_previous = 1.0e99;  // large value to force at least two rounds
    m_chi2 = 0.0;
    m_order = ord;
    m_fix_min = fix_min;
    m_fix_max = fix_max;
    m_max_it = std::abs(max_it);
    m_do_multilayer_rt_scale = do_multilayer_rt_scale;
    m_multilayer_rt_difference = std::make_unique<MultilayerRtDifference>(10000);
 
    m_tracker = std::make_unique<CurvedPatRec>();
 
    if (m_order == 0) {
        throw std::runtime_error(
            Form("File: %s, Line: %d\nRtCalibrationCurved::init - Order of the correction polynomial must be >0!", __FILE__, __LINE__));
    }
 
    m_t_length = 1000.0;
    m_t_mean = 500.0;
    // default values, correct values will be set when the input r-t
    // has been given
 
    m_U = std::vector<CLHEP::HepVector>(m_order);
    m_U_weighted = std::vector<CLHEP::HepVector>(m_order);
    m_A = CLHEP::HepSymMatrix(m_order, 0);
    m_b = CLHEP::HepVector(m_order, 0);
    m_alpha = CLHEP::HepVector(m_order, 0);
 
    // correction function
    if (func_type < 1 || func_type > 3) {
        throw std::runtime_error(
            Form("File: %s, Line: %d\nRtCalibrationCurved::init - Illegal correction function type!", __FILE__, __LINE__));
    }
    switch (func_type) {
        case 1: m_base_function = std::make_unique<LegendrePolynomial>(); break;
        case 2: m_base_function = std::make_unique<ChebyshevPolynomial>(); break;
        case 3:
            if (m_order < 2) {
                throw std::runtime_error(
                    Form("File: %s, Line: %d\nRtCalibrationCurved::init - Order must be >2 for polygons! It is set to %i by the user.",
                         __FILE__, __LINE__, m_order));
            }
            std::vector<double> x(m_order);
            double bin_width = 2.0 / static_cast<double>(m_order - 1);
            for (unsigned int k = 0; k < m_order; k++) { x[k] = -1 + k * bin_width; }
            m_base_function = std::make_unique<PolygonBase>(x);
            break;
    }
 
    // monotony of r(t) //
    m_force_monotony = false;
 
    // parabolic extrapolations //
    m_do_parabolic_extrapolation = do_parabolic_extrapolation;
 
    // smoothing //
    m_do_smoothing = do_smoothing;
 
    m_M_track = CLHEP::HepSymMatrix(3);
    m_M_track_inverse = CLHEP::HepSymMatrix(3);
 
    return;
}

◆ iteration()

int RtCalibrationCurved::iteration ( ) const

get the number of the current iteration

Definition at line 57 of file RtCalibrationCurved.cxx.

57{ return m_iteration; }

◆ name()

virtual std::string MuonCalib::IMdtCalibration::name ( ) const

inlinevirtualinherited

returns name (region) of instance

Definition at line 49 of file IMdtCalibration.h.

49{ return m_name; }

MuonCalib::IMdtCalibration::m_name

std::string m_name

Definition IMdtCalibration.h:52

◆ noParabolicExtrapolation()

void MuonCalib::RtCalibrationCurved::noParabolicExtrapolation ( )

no parabolic extrapolation is done

◆ noSmoothing()

void RtCalibrationCurved::noSmoothing ( )

do not smoothen the r-t relationship after convergence

Definition at line 100 of file RtCalibrationCurved.cxx.

100{ m_do_smoothing = false; }

◆ numberOfSegments()

int RtCalibrationCurved::numberOfSegments ( ) const

get the number of segments which were passed to the algorithm

Definition at line 53 of file RtCalibrationCurved.cxx.

53{ return m_nb_segments; }

◆ numberOfSegmentsUsed()

int RtCalibrationCurved::numberOfSegmentsUsed ( ) const

get the number of segments which are used in the autocalibration

Definition at line 55 of file RtCalibrationCurved.cxx.

55{ return m_nb_segments_used; }

◆ performParabolicExtrapolation()

std::shared_ptr< RtRelationLookUp > RtCalibrationCurved::performParabolicExtrapolation	(	const bool &	min,
		const bool &	max,
		const IRtRelation &	in_rt )

private

Definition at line 910 of file RtCalibrationCurved.cxx.

                                                                                                               {
    // VARIABLES //
    RtParabolicExtrapolation rt_extrapolator;           // r-t extrapolator
    std::shared_ptr<RtRelationLookUp> rt_low, rt_high;  // pointers to the r-t relationships after extrapolation
    std::vector<SamplePoint> add_fit_point;             // additional r-t points used if r(0) or r(t_max) is fixed.
    // EXTRAPOLATE TO LARGE RADII //
    if (max) {
        add_fit_point.clear();
        if (m_fix_max) { add_fit_point.push_back(SamplePoint(in_rt.tUpper(), in_rt.radius(in_rt.tUpper()), 1.0)); }
        if (in_rt.radius(in_rt.tUpper()) < 16.0) {
            rt_high = std::make_shared<RtRelationLookUp>(rt_extrapolator.getRtWithParabolicExtrapolation(
                in_rt, in_rt.radius(in_rt.tUpper()) - 3.0, in_rt.radius(in_rt.tUpper()) - 1.0, in_rt.radius(in_rt.tUpper()),
                add_fit_point));
        } else {
            rt_high =
                std::make_shared<RtRelationLookUp>(rt_extrapolator.getRtWithParabolicExtrapolation(in_rt, 13., 15., 16., add_fit_point));
        }
    }
    // EXTRAPOLATE TO SMALL RADII //
    if (min) {
        add_fit_point.clear();
        if (m_fix_min) { add_fit_point.push_back(SamplePoint(m_rt_new->tLower(), 0.0, 1.0)); }
        if (m_fix_max && rt_high != nullptr) {
            rt_low =
                std::make_shared<RtRelationLookUp>(rt_extrapolator.getRtWithParabolicExtrapolation(*rt_high, 1.0, 3.0, 0.0, add_fit_point));
        } else {
            rt_low =
                std::make_shared<RtRelationLookUp>(rt_extrapolator.getRtWithParabolicExtrapolation(in_rt, 1.0, 3.0, 0.0, add_fit_point));
        }
    }
    // RETURN RESULTS //
    if (min && max) {
        if (in_rt.hasTmaxDiff()) { rt_low->SetTmaxDiff(in_rt.GetTmaxDiff()); }
        return rt_low;
    }
    if (min) {
        if (in_rt.hasTmaxDiff()) { rt_low->SetTmaxDiff(in_rt.GetTmaxDiff()); }
        return rt_low;
    }
    if (in_rt.hasTmaxDiff() && rt_high) { rt_high->SetTmaxDiff(in_rt.GetTmaxDiff()); }
    return rt_high;
}

◆ reliability()

double RtCalibrationCurved::reliability ( ) const

get the reliability of the final r-t relationship: 0: no convergence yet 1: convergence, r(t) is reliable 2: convergence, r(t) is unreliable

Definition at line 49 of file RtCalibrationCurved.cxx.

49{ return m_status; }

◆ setEstimateRtAccuracy()

void RtCalibrationCurved::setEstimateRtAccuracy ( const double acc )

set the estimated r-t accuracy =acc

Definition at line 61 of file RtCalibrationCurved.cxx.

61{ m_rt_accuracy = std::abs(acc); }

◆ setInput()

void RtCalibrationCurved::setInput ( const IMdtCalibrationOutput * rt_input )

overridevirtual

set the r-t relationship, the internal autocalibration objects are reset

Implements MuonCalib::IMdtCalibration.

Definition at line 506 of file RtCalibrationCurved.cxx.

                                                                        {
    // VARIABLES //
    const RtCalibrationOutput *input(dynamic_cast<const RtCalibrationOutput *>(rt_input));
    // CHECK IF THE OUTPUT CLASS IS SUPPORTED //
    if (input == nullptr) {
        throw std::runtime_error(
            Form("File: %s, Line: %d\nRtCalibrationCurved::setInput - Calibration input class not supported.", __FILE__, __LINE__));
    }
    // GET THE INITIAL r-t RELATIONSHIP AND RESET STATUS VARIABLES //
 
    // get the r-t relationship //
    m_rt = input->rt();
 
    // status variables //
    m_nb_segments = 0;
    m_nb_segments_used = 0;
    m_chi2 = 0.0;
    m_A = CLHEP::HepSymMatrix(m_order, 0);
    m_b = CLHEP::HepVector(m_order, 0);
    m_alpha = CLHEP::HepVector(m_order, 0);
 
    // drift-time interval //
    auto rt_Chebyshev = dynamic_cast<const RtChebyshev*>(m_rt.get());
    auto rt_LookUp = dynamic_cast<const RtRelationLookUp*>(m_rt.get());
 
    if (!rt_Chebyshev && !rt_LookUp) {
        throw std::runtime_error(Form("File: %s, Line: %d\nRtCalibrationCurved::setInput - r-t class not supported.", __FILE__, __LINE__));
    }
 
    // RtChebyshev //
    if (rt_Chebyshev) {
        m_t_length = rt_Chebyshev->tUpper() - rt_Chebyshev->tLower();
        m_t_mean = 0.5 * (rt_Chebyshev->tLower() + rt_Chebyshev->tUpper());
    }
 
    // RtRelationLookUp, dangerous implementation, but the only way right now //
    if (rt_LookUp) {
        m_t_length = rt_LookUp->par(1) * (rt_LookUp->nPar() - 2) - rt_LookUp->par(0);
        m_t_mean = 0.5 * (rt_LookUp->par(1) * (rt_LookUp->nPar() - 2) + rt_LookUp->par(0));
    }
 
    return;
}

◆ smoothing()

bool RtCalibrationCurved::smoothing ( ) const

returns true, if the r-t relationship will be smoothened using the conventional autocalibration after convergence; returns false otherwise

Definition at line 59 of file RtCalibrationCurved.cxx.

59{ return m_do_smoothing; }

◆ switch_off_control_histograms()

void RtCalibrationCurved::switch_off_control_histograms ( )

the algorithm does not produce controll histograms (this is the default)

Definition at line 88 of file RtCalibrationCurved.cxx.

                                                        {
    m_control_histograms = false;
    if (m_tfile) {
        m_tfile->Write();
        m_multilayer_rt_difference = std::make_unique<MultilayerRtDifference>(10000);
    }
}

◆ switch_on_control_histograms()

void RtCalibrationCurved::switch_on_control_histograms ( const std::string & file_name )

this methods requests control histograms from the algorithms; the algorithm will write them to ROOT file called "file_name"

Definition at line 63 of file RtCalibrationCurved.cxx.

                                                                                 {
    // CREATE THE ROOT FILE AND THE HISTOGRAMS //
    m_control_histograms = true;
 
    m_tfile = std::make_unique<TFile>(file_name.c_str(), "RECREATE");
 
    m_cut_evolution = std::make_unique<TH1F>("m_cut_evolution", "CUT EVOLUTION", 11, -0.5, 10.5);
    m_nb_segment_hits = std::make_unique<TH1F>("m_nb_segment_hits", "NUMBER OF HITS ON THE REFITTED SEGMENTS", 11, -0.5, 10.5);
    m_pull_initial = std::make_unique<TH1F>("m_pull_initial", "INITIAL PULL DISTRIBUTION", 200, -5.05, 5.05);
    m_residuals_initial =
        std::make_unique<TH2F>("m_residuals_initial", "INITIAL OF THE REFITTED SEGMENTS", 100, -0.5, 15.0, 300, -1.5, 1.5);
    m_residuals_initial_all = std::make_unique<TH2F>("m_residuals_initial_all", "INITIAL OF THE REFITTED SEGMENTS BEFORE CONVERGENCE", 300,
                                                     -15.0, 15.0, 1000, -5, 5);
    m_residuals_final = std::make_unique<TH2F>("m_residuals_final", "FINAL OF THE REFITTED SEGMENTS", 100, -0.5, 15.0, 300, -1.5, 1.5);
    m_driftTime_final =
        std::make_unique<TH2F>("m_driftTime_final", "FINAL DRIFTTIME OF THE REFITTED SEGMENTS", 300, -15.0, 15.0, 300, -15.0, 15.0);
    m_driftTime_initial =
        std::make_unique<TH2F>("m_driftTime_initial", "FINAL DRIFTTIME OF THE REFITTED SEGMENTS", 300, -15.0, 15.0, 1300, -100, 1200);
    m_adc_vs_residual_final =
        std::make_unique<TH2F>("m_adc_resi_initial", "FINAL ADC VS RESIDUAL OF THE REFITTED SEGMENTS", 1350, 0, 1350, 300, -15, 15);
 
    m_multilayer_rt_difference = std::make_unique<MultilayerRtDifference>(10000, m_tfile.get());
}

◆ t_from_r()

double RtCalibrationCurved::t_from_r ( const double r )

private

Definition at line 803 of file RtCalibrationCurved.cxx.

                                                   {
    // VARIABLES //
    double precision(0.001);                            // spatial precision of the inversion
    double t_max(0.5 * (m_t_length + 2.0 * m_t_mean));  // upper time search limit
    double t_min(t_max - m_t_length);                   // lower time search limit
    // SEARCH FOR THE CORRESPONDING DRIFT TIME //
    while (t_max - t_min > 0.1 && std::abs(m_rt->radius(0.5 * (t_min + t_max)) - r) > precision) {
        if (m_rt->radius(0.5 * (t_min + t_max)) > r) {
            t_max = 0.5 * (t_min + t_max);
        } else {
            t_min = 0.5 * (t_min + t_max);
        }
    }
 
    return 0.5 * (t_min + t_max);
}

Member Data Documentation

◆ m_A

CLHEP::HepSymMatrix MuonCalib::RtCalibrationCurved::m_A

private

Definition at line 248 of file RtCalibrationCurved.h.

◆ m_adc_vs_residual_final

std::unique_ptr<TH2F> MuonCalib::RtCalibrationCurved::m_adc_vs_residual_final {}

private

Definition at line 268 of file RtCalibrationCurved.h.

268{}; // final residual distribution after convergence

◆ m_alpha

CLHEP::HepVector MuonCalib::RtCalibrationCurved::m_alpha

private

Definition at line 250 of file RtCalibrationCurved.h.

◆ m_b

CLHEP::HepVector MuonCalib::RtCalibrationCurved::m_b

private

Definition at line 252 of file RtCalibrationCurved.h.

◆ m_base_function

std::unique_ptr<BaseFunction> MuonCalib::RtCalibrationCurved::m_base_function

private

Definition at line 255 of file RtCalibrationCurved.h.

◆ m_chi2

double MuonCalib::RtCalibrationCurved::m_chi2 = 0.0

private

Definition at line 211 of file RtCalibrationCurved.h.

◆ m_chi2_previous

double MuonCalib::RtCalibrationCurved::m_chi2_previous = 0.0

private

Definition at line 205 of file RtCalibrationCurved.h.

◆ m_control_histograms

bool MuonCalib::RtCalibrationCurved::m_control_histograms = false

private

Definition at line 182 of file RtCalibrationCurved.h.

◆ m_cut_evolution

std::unique_ptr<TH1F> MuonCalib::RtCalibrationCurved::m_cut_evolution {}

private

Definition at line 259 of file RtCalibrationCurved.h.

259{}; // cut evolution histogram

◆ m_do_multilayer_rt_scale

bool MuonCalib::RtCalibrationCurved::m_do_multilayer_rt_scale = false

private

Definition at line 189 of file RtCalibrationCurved.h.

◆ m_do_parabolic_extrapolation

bool MuonCalib::RtCalibrationCurved::m_do_parabolic_extrapolation = false

private

Definition at line 234 of file RtCalibrationCurved.h.

◆ m_do_smoothing

bool MuonCalib::RtCalibrationCurved::m_do_smoothing = false

private

Definition at line 238 of file RtCalibrationCurved.h.

◆ m_driftTime_final

std::unique_ptr<TH2F> MuonCalib::RtCalibrationCurved::m_driftTime_final {}

private

Definition at line 267 of file RtCalibrationCurved.h.

267{}; // final residual distribution after convergence

◆ m_driftTime_initial

std::unique_ptr<TH2F> MuonCalib::RtCalibrationCurved::m_driftTime_initial {}

private

Definition at line 266 of file RtCalibrationCurved.h.

266{}; // final residual distribution after convergence

◆ m_fix_max

bool MuonCalib::RtCalibrationCurved::m_fix_max = false

private

Definition at line 185 of file RtCalibrationCurved.h.

◆ m_fix_min

bool MuonCalib::RtCalibrationCurved::m_fix_min = false

private

Definition at line 184 of file RtCalibrationCurved.h.

◆ m_force_monotony

bool MuonCalib::RtCalibrationCurved::m_force_monotony = false

private

Definition at line 187 of file RtCalibrationCurved.h.

◆ m_iteration

int MuonCalib::RtCalibrationCurved::m_iteration = 0

private

Definition at line 194 of file RtCalibrationCurved.h.

◆ m_M_track

CLHEP::HepSymMatrix MuonCalib::RtCalibrationCurved::m_M_track

private

Definition at line 230 of file RtCalibrationCurved.h.

◆ m_M_track_inverse

CLHEP::HepSymMatrix MuonCalib::RtCalibrationCurved::m_M_track_inverse

private

Definition at line 231 of file RtCalibrationCurved.h.

◆ m_max_it

int MuonCalib::RtCalibrationCurved::m_max_it = 0

private

Definition at line 186 of file RtCalibrationCurved.h.

◆ m_multilayer

std::array<bool, 2> MuonCalib::RtCalibrationCurved::m_multilayer {}

private

Definition at line 195 of file RtCalibrationCurved.h.

195{}; // m_multilayer[k] = true, if there was a segment

◆ m_multilayer_rt_difference

std::unique_ptr<MultilayerRtDifference> MuonCalib::RtCalibrationCurved::m_multilayer_rt_difference

private

Definition at line 224 of file RtCalibrationCurved.h.

◆ m_name

std::string MuonCalib::IMdtCalibration::m_name

privateinherited

Definition at line 52 of file IMdtCalibration.h.

◆ m_nb_segment_hits

std::unique_ptr<TH1F> MuonCalib::RtCalibrationCurved::m_nb_segment_hits {}

private

Definition at line 260 of file RtCalibrationCurved.h.

260{}; // number of hits on the segments

◆ m_nb_segments

int MuonCalib::RtCalibrationCurved::m_nb_segments = 0

private

Definition at line 192 of file RtCalibrationCurved.h.

◆ m_nb_segments_used

int MuonCalib::RtCalibrationCurved::m_nb_segments_used = 0

private

Definition at line 193 of file RtCalibrationCurved.h.

◆ m_order

unsigned int MuonCalib::RtCalibrationCurved::m_order = 0U

private

Definition at line 241 of file RtCalibrationCurved.h.

◆ m_output

std::shared_ptr<RtCalibrationOutput> MuonCalib::RtCalibrationCurved::m_output

private

Definition at line 222 of file RtCalibrationCurved.h.

◆ m_pull_final

std::unique_ptr<TH1F> MuonCalib::RtCalibrationCurved::m_pull_final {}

private

Definition at line 262 of file RtCalibrationCurved.h.

262{}; // final pull distribution after convergence

◆ m_pull_initial

std::unique_ptr<TH1F> MuonCalib::RtCalibrationCurved::m_pull_initial {}

private

Definition at line 261 of file RtCalibrationCurved.h.

261{}; // initial pull distribution

◆ m_r_max

double MuonCalib::RtCalibrationCurved::m_r_max = 0.0

private

Definition at line 226 of file RtCalibrationCurved.h.

◆ m_residuals_final

std::unique_ptr<TH2F> MuonCalib::RtCalibrationCurved::m_residuals_final {}

private

Definition at line 265 of file RtCalibrationCurved.h.

265{}; // final residual distribution after convergence

◆ m_residuals_initial

std::unique_ptr<TH2F> MuonCalib::RtCalibrationCurved::m_residuals_initial {}

private

Definition at line 263 of file RtCalibrationCurved.h.

263{}; // initial residual distribution

◆ m_residuals_initial_all

std::unique_ptr<TH2F> MuonCalib::RtCalibrationCurved::m_residuals_initial_all {}

private

Definition at line 264 of file RtCalibrationCurved.h.

264{}; // initial residual distribution before convergence

◆ m_rt

std::shared_ptr<const IRtRelation> MuonCalib::RtCalibrationCurved::m_rt

private

Definition at line 216 of file RtCalibrationCurved.h.

◆ m_rt_accuracy

double MuonCalib::RtCalibrationCurved::m_rt_accuracy = 0.0

private

Definition at line 202 of file RtCalibrationCurved.h.

◆ m_rt_accuracy_previous

double MuonCalib::RtCalibrationCurved::m_rt_accuracy_previous

private

Definition at line 203 of file RtCalibrationCurved.h.

◆ m_rt_new

std::shared_ptr<IRtRelation> MuonCalib::RtCalibrationCurved::m_rt_new

private

Definition at line 221 of file RtCalibrationCurved.h.

◆ m_status

int MuonCalib::RtCalibrationCurved::m_status = 0

private

Definition at line 199 of file RtCalibrationCurved.h.

◆ m_t_length

double MuonCalib::RtCalibrationCurved::m_t_length = 0.0

private

Definition at line 217 of file RtCalibrationCurved.h.

◆ m_t_mean

double MuonCalib::RtCalibrationCurved::m_t_mean = 0.0

private

Definition at line 218 of file RtCalibrationCurved.h.

◆ m_tfile

std::unique_ptr<TFile> MuonCalib::RtCalibrationCurved::m_tfile {}

private

Definition at line 258 of file RtCalibrationCurved.h.

258{}; // ROOT file

◆ m_tracker

std::unique_ptr<CurvedPatRec> MuonCalib::RtCalibrationCurved::m_tracker

private

Definition at line 227 of file RtCalibrationCurved.h.

◆ m_U

std::vector<CLHEP::HepVector> MuonCalib::RtCalibrationCurved::m_U

private

Definition at line 243 of file RtCalibrationCurved.h.

◆ m_U_weighted

std::vector<CLHEP::HepVector> MuonCalib::RtCalibrationCurved::m_U_weighted

private

Definition at line 244 of file RtCalibrationCurved.h.

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

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ MdtCalibOutputPtr

◆ MuonSegCit

◆ MuonSegIt

◆ MuonSegVec

Constructor & Destructor Documentation

◆ RtCalibrationCurved() [1/2]

◆ RtCalibrationCurved() [2/2]

◆ ~RtCalibrationCurved()

Member Function Documentation

◆ analyse()

◆ analyseSegments()

◆ converged()

◆ display_segment()

◆ doNotForceMonotony()

◆ doParabolicExtrapolation()

◆ doSmoothing()

◆ estimatedRtAccuracy()

◆ forceMonotony()

◆ fullMatrix()

◆ getResults()

◆ handleSegment()

◆ init()

◆ iteration()

◆ name()

◆ noParabolicExtrapolation()

◆ noSmoothing()

◆ numberOfSegments()

◆ numberOfSegmentsUsed()

◆ performParabolicExtrapolation()

◆ reliability()

◆ setEstimateRtAccuracy()

◆ setInput()

◆ smoothing()

◆ switch_off_control_histograms()

◆ switch_on_control_histograms()

◆ t_from_r()

Member Data Documentation

◆ m_A

◆ m_adc_vs_residual_final

◆ m_alpha

◆ m_b

◆ m_base_function

◆ m_chi2

◆ m_chi2_previous

◆ m_control_histograms

◆ m_cut_evolution

◆ m_do_multilayer_rt_scale

◆ m_do_parabolic_extrapolation

◆ m_do_smoothing

◆ m_driftTime_final

◆ m_driftTime_initial

◆ m_fix_max

◆ m_fix_min

◆ m_force_monotony

◆ m_iteration

◆ m_M_track

◆ m_M_track_inverse

◆ m_max_it

◆ m_multilayer

◆ m_multilayer_rt_difference

◆ m_name

◆ m_nb_segment_hits

◆ m_nb_segments

◆ m_nb_segments_used

◆ m_order

◆ m_output

◆ m_pull_final

◆ m_pull_initial

◆ m_r_max

◆ m_residuals_final

◆ m_residuals_initial

◆ m_residuals_initial_all

◆ m_rt

◆ m_rt_accuracy