This class performs the analytic autocalibration whose basic ideas were developed by Mario Deile (see ATL-MUON-2004-021). More...

#include <RtCalibrationAnalytic.h>

Inheritance diagram for MuonCalib::RtCalibrationAnalytic:

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
	RtCalibrationAnalytic (const std::string &name)
	Default constructor: r-t accuracy is set to 0.5 mm.
	RtCalibrationAnalytic (const std::string &name, const double rt_accuracy, const unsigned int &func_type, const unsigned int &ord, const bool &split, const bool &full_matrix, const bool &fix_min, const bool &fix_max, const int &max_it, bool do_smoothing=false, bool do_parabolic_extrapolation=false)
	Constructor: r-t accuracy is set to rt_accuracy (unit: CLHEP::mm).
	~RtCalibrationAnalytic ()
double	reliability () const
	get the reliability of the r-t: 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	splitIntoMultilayers () const
	returns true, if segments are internally restricted to single multilayers; returns false, if segments may contain hits from both multilayers
bool	fullMatrix () const
	returns true, if the full matrix relating the errors in the r-t relationship to the residuals should be used; returns false, if the unit matrix is used
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	splitIntoMultilayers (const bool &yes_or_no)
	yes_or_no=true: segments are internally restriced to single multilayers; yes_or_no=false: segments over two multilayers of a chamber are allowed
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
void	doNotForceMonotony ()
	do not force r(t) to be monotonically increasing
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
void	doParabolicExtrapolation ()
	requires that parabolic extrapolation will be used for small and large radii
void	noParabolicExtrapolation ()
	no parabolic extrapolation is done
MdtCalibOutputPtr	analyseSegments (const MuonSegVec &seg)
	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)
	set the r-t relationship, the internal autocalibration objects are reset
bool	analyse ()
	perform the autocalibration with the segments acquired so far
bool	converged () const
	returns true, if the autocalibration has converged
MdtCalibOutputPtr	getResults () const
	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 &split, const bool &full_matrix, const bool &fix_min, const bool &fix_max, const int &max_it, bool do_smoothing, bool do_parabolic_extrapolation)
double	t_from_r (const double r)
void	display_segment (MuonCalibSegment *segment, std::ofstream &outfile)
std::shared_ptr< RtRelationLookUp >	performParabolicExtrapolation (const bool &min, const bool &max, const IRtRelation &in_rt)

Private Attributes
bool	m_control_histograms = false
bool	m_split_into_ml = false
bool	m_full_matrix = false
bool	m_fix_min = false
bool	m_fix_max = false
int	m_max_it = 0
bool	m_force_monotony = 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 = 0.0
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
double	m_r_max = 0.0
QuasianalyticLineReconstruction	m_tracker
bool	m_do_smoothing = false
bool	m_do_parabolic_extrapolation = false
unsigned int	m_order = 0U
std::vector< CLHEP::HepVector >	m_U
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_CL
std::unique_ptr< TH2F >	m_residuals
MeanRMS	m_track_slope
MeanRMS	m_track_position
std::string	m_name

Detailed Description

This class performs the analytic autocalibration whose basic ideas were developed by Mario Deile (see ATL-MUON-2004-021).

Author: Olive.nosp@m.r.Ko.nosp@m.rtner.nosp@m.@CER.nosp@m.N.CH

Date: 06.04.2006

Author: Olive.nosp@m.r.Ko.nosp@m.rtner.nosp@m.@CER.nosp@m.N.CH

Definition at line 59 of file RtCalibrationAnalytic.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

◆ RtCalibrationAnalytic() [1/2]

RtCalibrationAnalytic::RtCalibrationAnalytic ( 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 5. Segments are restricted to single multilayers. The full matrix relating the errors in r(t) to the residuals is used. By default no smoothing is applied after convergence.

Definition at line 38 of file RtCalibrationAnalytic.cxx.

                                                                  : IMdtCalibration(name) {
    init(0.5 * CLHEP::mm, 1, 5, true, true, true, true, 100, false, false);
}

◆ RtCalibrationAnalytic() [2/2]

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

Constructor: r-t accuracy is set to rt_accuracy (unit: CLHEP::mm).

The r-t accuracy is used internally to distinguish between good and bad segments.

Parameters

func_type	type of function to be used for the r-t correction; = 1: Legendre polynomial, = 2: Chebyshev polynomial, = 3: polygon equidistant in r.
ord	The order of the r-t correction polynomial is set to ord.
split	= true forces the algorithm to restrict segments to multilayers. Otherwise, segments may contain hits from different multilayers.
full_matrix	= true lets the algorithm use the full matrix relating the errors in r(t) to the residuals is used. Otherwise the unit matrix is used, making the code equivalent to the conventional/classical method.
fix_min,fix_max	=true fix r(t_min), r(t_max) (this is default).
max_it	maximum number of iterations.
do_smoothing	Smoothen the r-t relations after convergence.
do_parabolic_extrapolation	Use parabolic extrapolations for small and larged drift radii.

Definition at line 42 of file RtCalibrationAnalytic.cxx.

                                                                                                                                         :
    IMdtCalibration(name), m_rt(nullptr) {
    init(rt_accuracy, func_type, ord, split, full_matrix, fix_min, fix_max, max_it, do_smoothing, do_parabolic_extrapolation);
}

◆ ~RtCalibrationAnalytic()

RtCalibrationAnalytic::~RtCalibrationAnalytic ( )

Definition at line 49 of file RtCalibrationAnalytic.cxx.

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

Member Function Documentation

◆ analyse()

bool RtCalibrationAnalytic::analyse ( )

perform the autocalibration with the segments acquired so far

Definition at line 740 of file RtCalibrationAnalytic.cxx.

                                    {
    if (m_tfile) m_tfile->cd();
    // 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
    RtParabolicExtrapolation rt_extrapolator;    // r-t extrapolator
    // SOLVE THE AUTOCALIBRATION EQUATION //
    m_alpha = m_A.inverse(ifail) * m_b;
    if (ifail != 0) {
        MsgStream log(Athena::getMessageSvc(), "RtCalibrationAnalytic");
        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 //
        nb_points = std::max(rt_Chebyshev->nDoF(), 30u);
 
        // 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() > 14.1) && 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()); }
            }
        }
 
        if (m_control_histograms) {
            std::unique_ptr<TGraphErrors> gre = std::make_unique<TGraphErrors>(x_r.size());
            for (unsigned int i = 0; i < 1; i++) {
                gre->SetPoint(i, x_r[i].x1(), x_r[i].x2());
                gre->SetPointError(i, 0, x_r[i].error());
            }
            std::ostringstream str_str;
            str_str << "CorrectionPoints_" << m_iteration;
            gre->Write(str_str.str().c_str());
        }
 
        // 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_unique<RtChebyshev>(rt_param);
    }
 
    // 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; }
 
        std::unique_ptr<TGraph> gr;
        if (m_control_histograms) gr = std::make_unique<TGraph>(max_k - min_k);
 
        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_control_histograms) gr->SetPoint(k - min_k, x, 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_unique<RtRelationLookUp>(rt_param);
        if (m_control_histograms) {
            std::ostringstream str_str;
            str_str << "CorrectionPoints_" << m_iteration;
            gr->Write(str_str.str().c_str());
        }
        m_r_max = m_rt_new->radius(m_rt_new->tUpper());
    }
    // DETERMINE THE r-t QUALITY AND CHECK FOR CONVERGENCE //
 
    // estimate r-t accuracy //
    m_rt_accuracy = 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_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_split_into_ml && m_nb_segments_used < 0.5 * (m_multilayer[0] + m_multilayer[1]) * m_nb_segments) { m_status = 2; }
        if (!m_split_into_ml && m_nb_segments_used < 0.5 * m_nb_segments) { m_status = 2; }
    }
    // STORE THE RESULTS IN THE RtCalibrationOutput //
    m_iteration = m_iteration + 1;
 
    m_output = std::make_shared<RtCalibrationOutput>(
        m_rt_new, std::make_shared<RtFullInfo>("RtCalibrationAnalytic", m_iteration, m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
 
    return true;
}

◆ analyseSegments()

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

virtual

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

Implements MuonCalib::IMdtCalibration.

Definition at line 339 of file RtCalibrationAnalytic.cxx.

                                                                                                   {
    std::shared_ptr<const IRtRelation> tmp_rt;
    std::shared_ptr<const IRtRelation> conv_rt;
    // AUTOCALIBRATION LOOP //
    while (!converged()) {
        for (const auto & k : seg) { handleSegment(*k); }
        if (!analyse()) {
            MsgStream log(Athena::getMessageSvc(), "RtCalibrationAnalytic");
            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();
 
        std::vector<double> params;
        params.push_back(tmp_rt->tLower());
        params.push_back(0.01 * (tmp_rt->tUpper() - tmp_rt->tLower()));
        for (double t = tmp_rt->tLower(); t <= tmp_rt->tUpper(); t = t + params[1]) { params.push_back(tmp_rt->radius(t)); }
        conv_rt = std::make_shared<RtRelationLookUp>(params);
    }
 
    // parabolic extrapolations for small radii //
    if (m_do_parabolic_extrapolation) {
        std::shared_ptr<const RtRelationLookUp> tmprt = performParabolicExtrapolation(true, true, *tmp_rt);
        m_output = std::make_shared<RtCalibrationOutput>(
            tmprt, std::make_shared<RtFullInfo>("RtCalibrationAnalyticExt", 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) { 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() < 3) { 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.1;
                }
            }
            avres = avres / static_cast<double>(k->mdtHitsOnTrack());
            avres = std::sqrt(avres);
            if (k->mdtHitsOnTrack() > 3) {
                if (smoothing.addResidualsFromSegment(*k, true, 7.0 * avres)) { counter++; }
            } else {
                if (smoothing.addResidualsFromSegment(*k, false, 7.0 * avres)) { counter++; }
            }
        }
 
        // break, do no smoothing if there are not enough segments //
        if (counter < 1000) {
            MsgStream log(Athena::getMessageSvc(), "RtCalibrationAnalytic");
            log << MSG::WARNING << "analyseSegments() - no smoothing applied due to too small number of reconstructed segments" << endmsg;
            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>("RtCalibrationAnalytic", m_iteration, m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
    // RETURN THE RESULT AFTER CONVERGENCE //
    return getResults();
}

◆ converged()

bool RtCalibrationAnalytic::converged ( ) const

returns true, if the autocalibration has converged

Definition at line 905 of file RtCalibrationAnalytic.cxx.

905{ return (m_status > 0); }

◆ display_segment()

void RtCalibrationAnalytic::display_segment	(	MuonCalibSegment *	segment,
		std::ofstream &	outfile )

private

Definition at line 170 of file RtCalibrationAnalytic.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 //
    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";
 
    // add a wait statement //
    outfile << "WAIT\n";
 
    return;
}

◆ doNotForceMonotony()

void RtCalibrationAnalytic::doNotForceMonotony ( )

do not force r(t) to be monotonically increasing

Definition at line 331 of file RtCalibrationAnalytic.cxx.

331{ m_force_monotony = false; }

◆ doParabolicExtrapolation()

void RtCalibrationAnalytic::doParabolicExtrapolation ( )

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

Definition at line 334 of file RtCalibrationAnalytic.cxx.

334{ m_do_parabolic_extrapolation = true; }

◆ doSmoothing()

void RtCalibrationAnalytic::doSmoothing ( )

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

Definition at line 332 of file RtCalibrationAnalytic.cxx.

332{ m_do_smoothing = true; }

◆ estimatedRtAccuracy()

double RtCalibrationAnalytic::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 257 of file RtCalibrationAnalytic.cxx.

257{ return m_rt_accuracy; }

◆ forceMonotony()

void RtCalibrationAnalytic::forceMonotony ( )

force r(t) to be monotonically increasing

Definition at line 330 of file RtCalibrationAnalytic.cxx.

330{ m_force_monotony = true; }

◆ fullMatrix() [1/2]

bool RtCalibrationAnalytic::fullMatrix ( ) const

returns true, if the full matrix relating the errors in the r-t relationship to the residuals should be used; returns false, if the unit matrix is used

Definition at line 292 of file RtCalibrationAnalytic.cxx.

292{ return m_full_matrix; }

MuonCalib::RtCalibrationAnalytic::m_full_matrix

bool m_full_matrix

Definition RtCalibrationAnalytic.h:223

◆ fullMatrix() [2/2]

void RtCalibrationAnalytic::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

Definition at line 338 of file RtCalibrationAnalytic.cxx.

338{ m_full_matrix = yes_or_no; }

◆ getResults()

RtCalibrationAnalytic::MdtCalibOutputPtr RtCalibrationAnalytic::getResults ( ) const

virtual

returns the final r-t relationship

Implements MuonCalib::IMdtCalibration.

Definition at line 906 of file RtCalibrationAnalytic.cxx.

906{ return m_output; }

◆ handleSegment()

bool RtCalibrationAnalytic::handleSegment ( MuonCalibSegment & seg )

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

Definition at line 468 of file RtCalibrationAnalytic.cxx.

                                                               {
    // RETURN, IF THE SEGMENT HAD LESS THAN 3 HITS //
    if (m_control_histograms) { m_cut_evolution->Fill(0.0, 1.0); }
 
    m_nb_segments = m_nb_segments + 1;
    if (seg.mdtHitsOnTrack() < 3) { 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[2];
    hit_selection[0] = IMdtSegmentFitter::HitSelection(seg.mdtHitsOnTrack());
    hit_selection[1] = IMdtSegmentFitter::HitSelection(seg.mdtHitsOnTrack());
    // hit selection vectors for refits in the first and second multilayer
    unsigned int nb_hits_in_ml[2];       // 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<double> residual_value;  // residuals
    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 //
    CLHEP::HepVector G;        // vector used in the calculation of the autocalibration matrix
    std::vector<double> zeta;  // vector used in the calculation of G
    // PREPARATION FOR THE SEGMENT REFIT //
 
    // debug display //
    //  display_segment(&seg, display);
 
    // 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 multilayer, if requested                    //
    nb_hits_in_ml[0] = 0;
    nb_hits_in_ml[1] = 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.5) {
            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);
        //      seg.mdtHOT()[k]->setDriftRadius(m_rt->radius(
        //                  seg.mdtHOT()[k]->driftTime()),
        //                  seg.mdtHOT()[k]->sigmaDriftRadius());
 
        // average resolution in the segment //
        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;
        }
 
        // hit selection //
        (hit_selection[0])[k] = 0;
 
        // 1st multilayer, if splitting is requested //
        if (m_split_into_ml) { (hit_selection[0])[k] = seg.mdtHOT()[k]->identify().mdtMultilayer() == 2; }
 
        // 2nd multilayer, if splitting is requested //
        if (m_split_into_ml) { (hit_selection[1])[k] = seg.mdtHOT()[k]->identify().mdtMultilayer() == 1; }
 
        // reject hits with bad radii or bad resolution //
        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[0])[k] = 1;
            (hit_selection[1])[k] = 1;
        }
 
        // counting of selected segments //
        nb_hits_in_ml[0] = nb_hits_in_ml[0] + (1 - (hit_selection[0])[k]);
        nb_hits_in_ml[1] = nb_hits_in_ml[1] + (1 - (hit_selection[1])[k]);
 
        // check for hits in both multilayers (needed if splitting is requested) //
        if (m_multilayer[0] == false && seg.mdtHOT()[k]->identify().mdtMultilayer() == 1) { m_multilayer[0] = true; }
        if (m_multilayer[1] == false && seg.mdtHOT()[k]->identify().mdtMultilayer() == 2) { m_multilayer[1] = true; }
    }
 
    // average resolution and chi^2 scale factor //
    av_res = std::sqrt(av_res / static_cast<double>(seg.mdtHitsOnTrack()));
    chi2_scale_factor = std::sqrt(av_res * av_res + m_rt_accuracy * m_rt_accuracy) / av_res;
    // FILL THE AUTOCALIBRATION MATRICES //
 
    // set the road width for the track reconstruction //
    m_tracker.setRoadWidth(7.0 * std::sqrt(av_res * av_res + m_rt_accuracy * m_rt_accuracy));
 
    // loop over the multilayers //
 
    //-----------------------------------------------------------------------------
    for (int k = 0; k < 1 + m_split_into_ml; k++) {
        //-----------------------------------------------------------------------------
 
        if (nb_hits_in_ml[k] < 3) { continue; }
 
        // refit the segments within the multilayers //
        MTStraightLine track;
 
        if (!m_tracker.fit(seg, hit_selection[k], track)) { continue; }
 
        if (std::isnan(track.a_x1()) || std::isnan(track.a_x2()) || std::isnan(track.b_x1()) || std::isnan(track.b_x2())) { continue; }
 
        // check the quality of the fit //
        if (track.chi2PerDegreesOfFreedom() > 5 * chi2_scale_factor) { continue; }
 
        // check whether there are at least three track hits //
        if (m_control_histograms) { m_nb_segment_hits->Fill(track.numberOfTrackHits(), 1.0); }
        if (track.numberOfTrackHits() < 3) { continue; }
 
        // reject tracks with silly parameters
        if (std::abs(track.a_x2()) > 8e8) continue;
 
        // for filling into data class
        if (m_iteration == 0) {
            m_track_slope.Insert(track.a_x2());
            m_track_position.Insert(track.b_x2());
        }
 
        // bookkeeping for convergence decision and reliability estimation //
        m_chi2 = m_chi2 + track.chi2PerDegreesOfFreedom();
        m_nb_segments_used = m_nb_segments_used + 1;
 
        // fill the autocalibration objects //
        // auxiliary variables //
        d_track = std::vector<double>(track.numberOfTrackHits());
        residual_value = std::vector<double>(track.numberOfTrackHits());
        w = std::vector<MTStraightLine>(track.numberOfTrackHits());
        G = CLHEP::HepVector(track.numberOfTrackHits());
        zeta = std::vector<double>(track.numberOfTrackHits());
 
        // base function values //
        for (unsigned int l = 0; l < m_order; l++) {
            m_U[l] = CLHEP::HepVector(track.numberOfTrackHits());
            for (unsigned int m = 0; m < track.numberOfTrackHits(); m++) {
                x = (track.trackHits()[m]->driftRadius() - 0.5 * m_r_max) / (0.5 * m_r_max);
                (m_U[l])[m] = m_base_function->value(l, x);
            }
        }
 
        // get the wire coordinates, track distances, and residuals //
        for (unsigned int l = 0; l < track.numberOfTrackHits(); l++) {
            w[l] =
                MTStraightLine(Amg::Vector3D(0.0, (track.trackHits()[l]->localPosition()).y(), (track.trackHits()[l]->localPosition()).z()),
                               xhat, null, null);
            d_track[l] = track.signDistFrom(w[l]);
            residual_value[l] = track.trackHits()[l]->driftRadius() - std::abs(d_track[l]);
            if (m_control_histograms) { m_residuals->Fill(std::abs(d_track[l]), residual_value[l], 1.0); }
        }
 
        // loop over all track hits //
        for (unsigned int l = 0; l < track.numberOfTrackHits(); l++) {
            // analytic calculation of G //
            for (unsigned int m = 0; m < track.numberOfTrackHits(); m++) {
                zeta[m] = std::sqrt(1.0 + std::pow(track.a_x2(), 2)) * (w[m].positionVector()).z() - track.a_x2() * d_track[m];
            }
            for (unsigned int m = 0; m < track.numberOfTrackHits(); m++) {
                double sum1(0.0), sum2(0.0), sum3(0.0), sum4(0.0);
                for (unsigned int ll = 0; ll < track.numberOfTrackHits(); ll++) {
                    sum1 = sum1 + (zeta[l] - zeta[ll]) * (zeta[ll] - zeta[m]) /
                                      (track.trackHits()[m]->sigma2DriftRadius() * track.trackHits()[ll]->sigma2DriftRadius());
                    sum2 = sum2 + zeta[ll] / track.trackHits()[ll]->sigma2DriftRadius();
                    sum3 = sum3 + 1.0 / track.trackHits()[ll]->sigma2DriftRadius();
                    sum4 = sum4 + std::pow(zeta[ll] / track.trackHits()[ll]->sigmaDriftRadius(), 2);
                }
                if (d_track[m] * d_track[l] >= 0.0) {
                    G[m] = (l == m) - m_full_matrix * sum1 / (sum2 * sum2 - sum3 * sum4);
                } else {
                    G[m] = (l == m) + m_full_matrix * sum1 / (sum2 * sum2 - sum3 * sum4);
                }
            }
            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(G, m_U[p]) * dot(G, m_U[pp])) / track.trackHits()[l]->sigma2DriftRadius();
                    if (std::isnan(A_tmp[p][pp])) return true;
                }
            }
            m_A = A_tmp;
            CLHEP::HepVector b_tmp(m_b);
            for (unsigned int p = 0; p < m_order; p++) {
                b_tmp[p] = m_b[p] - residual_value[l] * dot(G, m_U[p]) / track.trackHits()[l]->sigma2DriftRadius();
                if (std::isnan(b_tmp[p])) return true;
            }
            m_b = b_tmp;
        }
 
        //-----------------------------------------------------------------------------
    }
    //-----------------------------------------------------------------------------
    return true;
}

◆ init()

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

private

Definition at line 56 of file RtCalibrationAnalytic.cxx.

                                                                  {
    // RESET PRIVATE VARIABLES //
    m_r_max = 15.0 * CLHEP::mm;
    m_control_histograms = false;
    m_tfile = nullptr;
    m_cut_evolution = nullptr;
    m_nb_segment_hits = nullptr;
    m_CL = nullptr;
    m_residuals = nullptr;
    m_split_into_ml = split;
    m_full_matrix = full_matrix;
    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_rt_accuracy_previous = 0.;
    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);
 
    if (m_order == 0) {
        throw std::runtime_error(
            Form("File: %s, Line: %d\nRtCalibrationAnalytic::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_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) {
        m_base_function = nullptr;
        throw std::runtime_error(
            Form("File: %s, Line: %d\nRtCalibrationAnalytic::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\nRtCalibrationAnalytic::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;
    }
 
    // request a chi^2 refit after the quasianalytic pattern recognition //
    m_tracker.switchOnRefit();
 
    // monotony of r(t) //
    m_force_monotony = false;
 
    // smoothing //
    m_do_smoothing = do_smoothing;
 
    // parabolic extrapolation //
    m_do_parabolic_extrapolation = do_parabolic_extrapolation;
 
    return;
}

◆ iteration()

int RtCalibrationAnalytic::iteration ( ) const

get the number of the current iteration

Definition at line 278 of file RtCalibrationAnalytic.cxx.

278{ 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 RtCalibrationAnalytic::noParabolicExtrapolation ( )

no parabolic extrapolation is done

Definition at line 335 of file RtCalibrationAnalytic.cxx.

335{ m_do_parabolic_extrapolation = false; }

◆ noSmoothing()

void RtCalibrationAnalytic::noSmoothing ( )

do not smoothen the r-t relationship after convergence

Definition at line 333 of file RtCalibrationAnalytic.cxx.

333{ m_do_smoothing = false; }

◆ numberOfSegments()

int RtCalibrationAnalytic::numberOfSegments ( ) const

get the number of segments which were passed to the algorithm

Definition at line 264 of file RtCalibrationAnalytic.cxx.

264{ return m_nb_segments; }

◆ numberOfSegmentsUsed()

int RtCalibrationAnalytic::numberOfSegmentsUsed ( ) const

get the number of segments which are used in the autocalibration

Definition at line 271 of file RtCalibrationAnalytic.cxx.

271{ return m_nb_segments_used; }

◆ performParabolicExtrapolation()

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

private

Definition at line 907 of file RtCalibrationAnalytic.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, m_r_max - 3.0, m_r_max - 1.0, m_r_max, 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) {
            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) { return rt_low; }
    if (min) { return rt_low; }
    return rt_high;
}

◆ reliability()

double RtCalibrationAnalytic::reliability ( ) const

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

Definition at line 250 of file RtCalibrationAnalytic.cxx.

250{ return m_status; }

◆ setEstimateRtAccuracy()

void RtCalibrationAnalytic::setEstimateRtAccuracy ( const double acc )

set the estimated r-t accuracy =acc

Definition at line 336 of file RtCalibrationAnalytic.cxx.

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

◆ setInput()

void RtCalibrationAnalytic::setInput ( const IMdtCalibrationOutput * rt_input )

virtual

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

Implements MuonCalib::IMdtCalibration.

Definition at line 688 of file RtCalibrationAnalytic.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\nRtCalibrationAnalytic::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();
    m_r_max = m_rt->radius(m_rt->tUpper());
 
    // 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\nRtCalibrationAnalytic::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));
    }
}

◆ smoothing()

bool RtCalibrationAnalytic::smoothing ( ) const

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

Definition at line 299 of file RtCalibrationAnalytic.cxx.

299{ return m_do_smoothing; }

◆ splitIntoMultilayers() [1/2]

bool RtCalibrationAnalytic::splitIntoMultilayers ( ) const

returns true, if segments are internally restricted to single multilayers; returns false, if segments may contain hits from both multilayers

Definition at line 285 of file RtCalibrationAnalytic.cxx.

285{ return m_split_into_ml; }

◆ splitIntoMultilayers() [2/2]

void RtCalibrationAnalytic::splitIntoMultilayers ( const bool & yes_or_no )

yes_or_no=true: segments are internally restriced to single multilayers; yes_or_no=false: segments over two multilayers of a chamber are allowed

Definition at line 337 of file RtCalibrationAnalytic.cxx.

337{ m_split_into_ml = yes_or_no; }

◆ switch_off_control_histograms()

void RtCalibrationAnalytic::switch_off_control_histograms ( )

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

Definition at line 322 of file RtCalibrationAnalytic.cxx.

                                                          {
    m_control_histograms = false;
    if (m_tfile) {
        m_tfile->Write();
        m_tfile.reset();
    }
}

◆ switch_on_control_histograms()

void RtCalibrationAnalytic::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 306 of file RtCalibrationAnalytic.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_CL = std::make_unique<TH1F>("m_CL", "CONFIDENCE LEVELS OF THE REFITTED SEGMENTS", 100, 0.0, 1.0);
 
    m_residuals = std::make_unique<TH2F>("m_residuals", "RESIDUALS OF THE REFITTED SEGMENTS", 100, -0.5, 15.0, 300, -1.5, 1.5);
}

◆ t_from_r()

double RtCalibrationAnalytic::t_from_r ( const double r )

private

Definition at line 143 of file RtCalibrationAnalytic.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::RtCalibrationAnalytic::m_A

private

Definition at line 288 of file RtCalibrationAnalytic.h.

◆ m_alpha

CLHEP::HepVector MuonCalib::RtCalibrationAnalytic::m_alpha

private

Definition at line 290 of file RtCalibrationAnalytic.h.

◆ m_b

CLHEP::HepVector MuonCalib::RtCalibrationAnalytic::m_b

private

Definition at line 292 of file RtCalibrationAnalytic.h.

◆ m_base_function

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

private

Definition at line 295 of file RtCalibrationAnalytic.h.

◆ m_chi2

double MuonCalib::RtCalibrationAnalytic::m_chi2 = 0.0

private

Definition at line 255 of file RtCalibrationAnalytic.h.

◆ m_chi2_previous

double MuonCalib::RtCalibrationAnalytic::m_chi2_previous = 0.0

private

Definition at line 249 of file RtCalibrationAnalytic.h.

◆ m_CL

std::unique_ptr<TH1F> MuonCalib::RtCalibrationAnalytic::m_CL

private

Definition at line 301 of file RtCalibrationAnalytic.h.

◆ m_control_histograms

bool MuonCalib::RtCalibrationAnalytic::m_control_histograms = false

private

Definition at line 217 of file RtCalibrationAnalytic.h.

◆ m_cut_evolution

std::unique_ptr<TH1F> MuonCalib::RtCalibrationAnalytic::m_cut_evolution

private

Definition at line 299 of file RtCalibrationAnalytic.h.

◆ m_do_parabolic_extrapolation

bool MuonCalib::RtCalibrationAnalytic::m_do_parabolic_extrapolation = false

private

Definition at line 281 of file RtCalibrationAnalytic.h.

◆ m_do_smoothing

bool MuonCalib::RtCalibrationAnalytic::m_do_smoothing = false

private

Definition at line 278 of file RtCalibrationAnalytic.h.

◆ m_fix_max

bool MuonCalib::RtCalibrationAnalytic::m_fix_max = false

private

Definition at line 230 of file RtCalibrationAnalytic.h.

◆ m_fix_min

bool MuonCalib::RtCalibrationAnalytic::m_fix_min = false

private

Definition at line 229 of file RtCalibrationAnalytic.h.

◆ m_force_monotony

bool MuonCalib::RtCalibrationAnalytic::m_force_monotony = false

private

Definition at line 232 of file RtCalibrationAnalytic.h.

◆ m_full_matrix

bool MuonCalib::RtCalibrationAnalytic::m_full_matrix = false

private

Definition at line 223 of file RtCalibrationAnalytic.h.

◆ m_iteration

int MuonCalib::RtCalibrationAnalytic::m_iteration = 0

private

Definition at line 238 of file RtCalibrationAnalytic.h.

◆ m_max_it

int MuonCalib::RtCalibrationAnalytic::m_max_it = 0

private

Definition at line 231 of file RtCalibrationAnalytic.h.

◆ m_multilayer

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

private

Definition at line 239 of file RtCalibrationAnalytic.h.

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

◆ 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::RtCalibrationAnalytic::m_nb_segment_hits

private

Definition at line 300 of file RtCalibrationAnalytic.h.

◆ m_nb_segments

int MuonCalib::RtCalibrationAnalytic::m_nb_segments = 0

private

Definition at line 236 of file RtCalibrationAnalytic.h.

◆ m_nb_segments_used

int MuonCalib::RtCalibrationAnalytic::m_nb_segments_used = 0

private

Definition at line 237 of file RtCalibrationAnalytic.h.

◆ m_order

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

private

Definition at line 285 of file RtCalibrationAnalytic.h.

◆ m_output

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

private

Definition at line 266 of file RtCalibrationAnalytic.h.

◆ m_r_max

double MuonCalib::RtCalibrationAnalytic::m_r_max = 0.0

private

Definition at line 270 of file RtCalibrationAnalytic.h.

◆ m_residuals

std::unique_ptr<TH2F> MuonCalib::RtCalibrationAnalytic::m_residuals

private

Definition at line 302 of file RtCalibrationAnalytic.h.

◆ m_rt

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

private

Definition at line 260 of file RtCalibrationAnalytic.h.

◆ m_rt_accuracy

double MuonCalib::RtCalibrationAnalytic::m_rt_accuracy = 0.0

private

Definition at line 246 of file RtCalibrationAnalytic.h.

◆ m_rt_accuracy_previous

double MuonCalib::RtCalibrationAnalytic::m_rt_accuracy_previous = 0.0

private

Definition at line 247 of file RtCalibrationAnalytic.h.

◆ m_rt_new

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

private

Definition at line 265 of file RtCalibrationAnalytic.h.

◆ m_split_into_ml

bool MuonCalib::RtCalibrationAnalytic::m_split_into_ml = false

private

Definition at line 219 of file RtCalibrationAnalytic.h.

◆ m_status

int MuonCalib::RtCalibrationAnalytic::m_status = 0

private

Definition at line 243 of file RtCalibrationAnalytic.h.

◆ m_t_length

double MuonCalib::RtCalibrationAnalytic::m_t_length = 0.0

private

Definition at line 261 of file RtCalibrationAnalytic.h.

◆ m_t_mean

double MuonCalib::RtCalibrationAnalytic::m_t_mean = 0.0

private

Definition at line 262 of file RtCalibrationAnalytic.h.

◆ m_tfile

std::unique_ptr<TFile> MuonCalib::RtCalibrationAnalytic::m_tfile

private

Definition at line 298 of file RtCalibrationAnalytic.h.

◆ m_track_position

MeanRMS MuonCalib::RtCalibrationAnalytic::m_track_position

private

Definition at line 346 of file RtCalibrationAnalytic.h.

◆ m_track_slope

MeanRMS MuonCalib::RtCalibrationAnalytic::m_track_slope

private

Definition at line 345 of file RtCalibrationAnalytic.h.

◆ m_tracker

QuasianalyticLineReconstruction MuonCalib::RtCalibrationAnalytic::m_tracker

private

Definition at line 271 of file RtCalibrationAnalytic.h.

◆ m_U

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

private

Definition at line 287 of file RtCalibrationAnalytic.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

◆ RtCalibrationAnalytic() [1/2]

◆ RtCalibrationAnalytic() [2/2]

◆ ~RtCalibrationAnalytic()

Member Function Documentation

◆ analyse()

◆ analyseSegments()

◆ converged()

◆ display_segment()

◆ doNotForceMonotony()

◆ doParabolicExtrapolation()

◆ doSmoothing()

◆ estimatedRtAccuracy()

◆ forceMonotony()

◆ fullMatrix() [1/2]

◆ fullMatrix() [2/2]

◆ getResults()

◆ handleSegment()

◆ init()

◆ iteration()

◆ name()

◆ noParabolicExtrapolation()

◆ noSmoothing()

◆ numberOfSegments()

◆ numberOfSegmentsUsed()

◆ performParabolicExtrapolation()

◆ reliability()

◆ setEstimateRtAccuracy()

◆ setInput()

◆ smoothing()

◆ splitIntoMultilayers() [1/2]

◆ splitIntoMultilayers() [2/2]

◆ switch_off_control_histograms()

◆ switch_on_control_histograms()

◆ t_from_r()

Member Data Documentation

◆ m_A

◆ m_alpha

◆ m_b

◆ m_base_function

◆ m_chi2

◆ m_chi2_previous

◆ m_CL

◆ m_control_histograms

◆ m_cut_evolution

◆ m_do_parabolic_extrapolation

◆ m_do_smoothing

◆ m_fix_max

◆ m_fix_min

◆ m_force_monotony

◆ m_full_matrix

◆ m_iteration

◆ m_max_it

◆ m_multilayer

◆ m_name

◆ m_nb_segment_hits

◆ m_nb_segments

◆ m_nb_segments_used

◆ m_order

◆ m_output

◆ m_r_max

◆ m_residuals

◆ m_rt

◆ m_rt_accuracy

◆ m_rt_accuracy_previous

◆ m_rt_new

◆ m_split_into_ml

◆ m_status

◆ m_t_length

◆ m_t_mean