RtCalibrationCurved.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef MuonCalib_RtCalibrationCurvedH
#define MuonCalib_RtCalibrationCurvedH

 
#include <list>
#include <memory>
#include <string>
#include <vector>
 
// CLHEP //
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Vector.h"
 
// ROOT //
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
 
// MuonCalib //
#include "GeoPrimitives/GeoPrimitives.h"
#include "MdtCalibInterfaces/IMdtCalibration.h"
 
namespace MuonCalib {
 
    class IMdtCalibrationOutput;
    class IRtRelation;
    class RtRelationLookUp;
    class RtCalibrationOutput;
    class CurvedPatRec;
    class MuonCalibSegment;
    class BaseFunction;
    class CurvedLine;
    class MultilayerRtDifference;
 
    class RtCalibrationCurved : public IMdtCalibration {
    public:
        // Constructors //
        RtCalibrationCurved(const std::string &name);
 
        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);
 
        // Destructor //
        ~RtCalibrationCurved();
 
        // Methods //
        // get-methods //
        double reliability() const;
        double estimatedRtAccuracy() const;
        int numberOfSegments() const;
        int numberOfSegmentsUsed() const;
        int iteration() const;
        bool smoothing() const;
 
        // set-method //
        void setEstimateRtAccuracy(const double acc);
        void fullMatrix(const bool &yes_or_no);
        void switch_on_control_histograms(const std::string &file_name);
        void switch_off_control_histograms();
        void forceMonotony();
        void doNotForceMonotony();
        void doParabolicExtrapolation();
        void noParabolicExtrapolation();
        void doSmoothing();
        void noSmoothing();
 
        // methods required by the base class "IMdtCalibration" //
        MdtCalibOutputPtr analyseSegments(const MuonSegVec &seg) override;
        bool handleSegment(MuonCalibSegment &seg);
        void setInput(const IMdtCalibrationOutput *rt_input) override;
        bool analyse(const MuonSegVec &seg);
        bool converged() const;
        virtual MdtCalibOutputPtr getResults() const override;
 
    private:
        // options //
        bool m_control_histograms = false;      // = true, if control histograms should be
                                        //         produces
        bool m_fix_min = false;                 // = true: fix r(t_min)
        bool m_fix_max = false;                 // = true: fix r(t_max)
        int m_max_it = 0;                   // maximum number of iterations
        bool m_force_monotony = false;          // = true if r(t) is forced to monotonically
                                        //   increasing, false otherwise
        bool m_do_multilayer_rt_scale = false;  // determine multilayer rt scaling
 
        // bookkeeping //
        int m_nb_segments = 0;                 // number of segments passed to the algorithm
        int m_nb_segments_used = 0;            // number of segments used by the algorithm
        int m_iteration = 0;                   // current iteration
        std::array<bool, 2> m_multilayer{};  // m_multilayer[k] = true, if there was a segment
                                           //                   extending to multilayer k+1
 
        // r-t quality //
        int m_status = 0;                   // m_status: 0: no covergence yet,
                                        //           1: convergence, r-t is reliable,
                                        //           2: convergence, r-t is unreliable
        double m_rt_accuracy = 0.0;           // r-t accuracy (CLHEP::mm) of the input r-t
        double m_rt_accuracy_previous;  // r-t accuracy of the previous iteration
                                        // (used in the convergence criterion)
        double m_chi2_previous = 0.0;
        // average chi^2 per degrees of freedom from the
        // previous iteration (set to a large initial value
        // to force at least two iterations);
        // if an iteration gives a larger average than the
        // pervious iteration, the algorithm has converged
        double m_chi2 = 0.0;  // average chi^2 per degrees of freedom,
                        // if an iteration gives a larger average than the
                        // pervious iteration, the algorithm has converged
 
        // r-t relationship //
        std::shared_ptr<const IRtRelation> m_rt;  // pointer to the input r-t relationship
        double m_t_length = 0.0;                        // size of the drift time interval
        double m_t_mean = 0.0;                          // mean value of the drift time interval
 
        // r-t output //
        std::shared_ptr<IRtRelation> m_rt_new;          // r-t as determined by the autocalibration
        std::shared_ptr<RtCalibrationOutput> m_output;  // class holding the results of the
                                                        // autocalibration
        std::unique_ptr<MultilayerRtDifference> m_multilayer_rt_difference;
        // curved-segment fitting //
        double m_r_max = 0.0;                           // maximum value for accepted drift radii
        std::unique_ptr<CurvedPatRec> m_tracker;  // curved segment finder (used for track fitting)
                                                  // The following three objects are needed for autocalibration formulae.
 
        CLHEP::HepSymMatrix m_M_track;          // segment parameters = m_M_track^-1*m_b_track
        CLHEP::HepSymMatrix m_M_track_inverse;  // inverse of m_M_track
 
        // autocalibration objects //
        bool m_do_parabolic_extrapolation = false;           // = true: parabolic extrapolation is
                                                     //         done for small and large radii
                                                     // = false: no parabolic extrapolation is
                                                     //          done
        bool m_do_smoothing = false;                         // = true: the r-t relationship is smoothened after
                                                     //         convergence, no smoothing is done
                                                     //         otherwise
        unsigned int m_order = 0U;                        // order of the polynomial describing the
                                                     // correction to the r-t relationship
        std::vector<CLHEP::HepVector> m_U;           // vector of base correction function values
        std::vector<CLHEP::HepVector> m_U_weighted;  // vector of base correction function
                                                     // values weighted by the inverse
                                                     // standard deviation of the radius
                                                     // measurements
        CLHEP::HepSymMatrix m_A;                     // coefficient matrix of the final autocalibration
                                                     // equation
        CLHEP::HepVector m_alpha;                    // vector of fit parameters, i.e. the coefficients
                                                     // of the correction polynomial
        CLHEP::HepVector m_b;                        // m_A*m_alpha = m_b (final autocalibration equation)
 
        // correction functions //
        std::unique_ptr<BaseFunction> m_base_function;  // pointer to the base function u
 
        // control histograms //
        std::unique_ptr<TFile> m_tfile{};                 // ROOT file
        std::unique_ptr<TH1F> m_cut_evolution{};          // cut evolution histogram
        std::unique_ptr<TH1F> m_nb_segment_hits{};        // number of hits on the segments
        std::unique_ptr<TH1F> m_pull_initial{};           // initial pull distribution
        std::unique_ptr<TH1F> m_pull_final{};             // final pull distribution after convergence
        std::unique_ptr<TH2F> m_residuals_initial{};      // initial residual distribution
        std::unique_ptr<TH2F> m_residuals_initial_all{};  // initial residual distribution before convergence
        std::unique_ptr<TH2F> m_residuals_final{};        // final residual distribution after convergence
        std::unique_ptr<TH2F> m_driftTime_initial{};      // final residual distribution after convergence
        std::unique_ptr<TH2F> m_driftTime_final{};        // final residual distribution after convergence
        std::unique_ptr<TH2F> m_adc_vs_residual_final{};  // final residual distribution after convergence
 
        // private methods //
        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);
        // initialization method:
        // rt_accuracy = estimated r-t accuracy,
        // func_type: type of function to be used for
        //            the r-t correction;
        //     = 1: Legendre polynomial,
        //     = 2: Chebyshev polynomial,
        //     = 3: polygon
        // ord = "order" ot the r-t correction function
        // split = true forces the algorithm to restrict
        // segments to multilayers;
        // fix_min, fix_max=true: fix r(t_min), r(t_max)
        // max_it: maximum number of iterations
        // do_parabolic_extrapolation: do or do not use
        //    parabolic extrapolations for small and large
        //    radii;
        // do_smoothing: smoothen the r-t relationship
        //    after convergence if do_smoothing = true
        double t_from_r(const double r);
        // get t(r) for the input r-t relationship,
        // the method is auxiliary and not optimized;
        // it will disappear when the t(r) will be
        // available in the MuonCalib framework
        void display_segment(MuonCalibSegment *segment, std::ofstream &outfile, const CurvedLine *curved_segment);
        // write out a simple PAW macro displaying the
        // segment; if the pointer "curved_segment" equals
        // 0, a straight line as stored in segment is drawn;
        // the curved line is used otherwise
        std::shared_ptr<RtRelationLookUp> performParabolicExtrapolation(const bool &min, const bool &max, const IRtRelation &in_rt);
        // use parabolic extrapolations on the given r-t
        // relationship in_rt;
        // min: if true, use parabolic extrapolation towards
        //      r=0;
        // max: if true, use parabolic extrapolation towards
        //      r=r_max;
    };
 
}  // namespace MuonCalib
 
#endif