From 794cfaf225db8e53cfe013b92b5b7cfae063eaf7 Mon Sep 17 00:00:00 2001
From: Jochen Meyer <Jochen.Meyer@cern.ch>
Date: Mon, 15 Sep 2014 16:54:45 +0200
Subject: [PATCH] fixing compiler warning (MdtCalibRt-02-00-01)
---
.../MdtCalibRt/AdaptiveResidualSmoothing.h | 103 +
.../MdtCalib/MdtCalibRt/MdtCalibRt/MeanRMS.h | 43 +
.../MdtCalibRt/RtCalibrationAnalytic.h | 398 ++++
.../MdtCalibRt/RtCalibrationCurved.h | 352 ++++
.../MdtCalibRt/RtCalibrationIntegration.h | 127 ++
.../MdtCalibRt/RtCalibrationOutput.h | 36 +
.../MdtCalibRt/RtParabolicExtrapolation.h | 86 +
.../MdtCalib/MdtCalibRt/cmt/requirements | 31 +
.../MdtCalib/MdtCalibRt/doc/mainpage.h | 40 +
.../src/AdaptiveResidualSmoothing.cxx | 484 +++++
.../MdtCalibRt/src/MultilayerRtDifference.cxx | 220 +++
.../MdtCalibRt/src/MultilayerRtDifference.h | 47 +
.../MdtCalibRt/src/RtCalibrationAnalytic.cxx | 1641 ++++++++++++++++
.../MdtCalibRt/src/RtCalibrationCurved.cxx | 1699 +++++++++++++++++
.../src/RtCalibrationIntegration.cxx | 416 ++++
.../src/RtParabolicExtrapolation.cxx | 248 +++
16 files changed, 5971 insertions(+)
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/AdaptiveResidualSmoothing.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/MeanRMS.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationAnalytic.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationIntegration.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationOutput.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtParabolicExtrapolation.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/cmt/requirements
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/doc/mainpage.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/AdaptiveResidualSmoothing.cxx
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.cxx
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.h
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationAnalytic.cxx
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationIntegration.cxx
create mode 100644 MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtParabolicExtrapolation.cxx
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/AdaptiveResidualSmoothing.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/AdaptiveResidualSmoothing.h
new file mode 100644
index 0000000000000..23d4ea0abc2a0
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/AdaptiveResidualSmoothing.h
@@ -0,0 +1,103 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 29.07.2008, AUTHOR: OLIVER KORTNER
+// Modified: 20.03.2009 by O. Kortner, additional method ,,performSmoothing''
+// with better performance than the old method added;
+// softer chi^2 cuts in addResidualsFromSegment.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef MuonCalib_AdaptiveResidualSmoothingH
+#define MuonCalib_AdaptiveResidualSmoothingH
+
+//:::::::::::::::::::::::::::::::::::::
+//:: CLASS AdaptiveResidualSmoothing ::
+//:::::::::::::::::::::::::::::::::::::
+
+/// \class AdaptiveResidualSmoothing
+///
+/// This class is used to modify the r-t relationship to give smooth residuals
+/// centred at 0. The binning of the residual distribution is adapted to the
+/// statistics. The user can set the number of entries per bin. The user can
+/// decide whether the track fits should be performed with a straight or
+/// parabolic curve.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 28.07.2008
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// STL //
+#include <vector>
+
+// MuonCalib //
+#include "MdtCalibData/RtRelationLookUp.h"
+#include "MdtCalibFitters/StraightPatRec.h"
+#include "MdtCalibFitters/CurvedPatRec.h"
+#include "MuonCalibMath/DataPoint.h"
+
+namespace MuonCalib {
+
+class IRtRelation;
+
+class AdaptiveResidualSmoothing {
+
+public:
+// Constructor //
+ AdaptiveResidualSmoothing(void);
+ ///< Default constructor.
+
+// Methods //
+ void clear(void);
+ ///< clear the memory of the class
+ void addResidual(const double & radius, const double & residual);
+ ///< add the residual at the given radius
+ bool addResidualsFromSegment(MuonCalibSegment & seg, bool curved,
+ double road_width);
+ ///< reconstruct the given segment and
+ ///< store the residuals; if curved is true
+ ///< a curved segment fit is performed,
+ ///< otherwise a straight segment is fitted
+ ///< to the drift radii; the user must
+ ///< set the road width used in the pattern
+ ///< recognition;
+ ///< returns true in case of success
+ RtRelationLookUp performSmoothing(const IRtRelation & rt_rel,
+ unsigned int nb_entries_per_bin,
+ bool fix_t_min, bool fix_t_max);
+ ///< use the stored residuals to improve
+ ///< the given r-t relationship to give
+ ///< smoother residuals; the user has to
+ ///< set the number of entries per radial
+ ///< bin; the user can request that the
+ ///< radii for the minimum and maximum
+ ///< drift time are untouched
+ RtRelationLookUp performSmoothing(const IRtRelation & rt_rel,
+ const bool & fix_t_min,
+ const bool & fix_t_max);
+ ///< use the stored residuals to improve
+ ///< the given r-t relationship to give
+ ///< smoother residuals; the user can
+ ///< request that the radii for the minimum
+ ///< and maximum drift time are untouched
+
+private:
+ std::vector<DataPoint> m_residual_point; // vector of residual points
+ StraightPatRec m_sfitter; // straight-line fitter
+ CurvedPatRec m_cfitter; // curved-line fitter
+ double t_from_r(const IRtRelation & rt_rel, const double & r);
+ // get t(r) for the given r-t relationship,
+ // the method is auxiliary and not optimized;
+ // it will disappear when the t(r) will be
+ // available in the MuonCalib framework
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/MeanRMS.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/MeanRMS.h
new file mode 100644
index 0000000000000..2f76af94d6fbe
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/MeanRMS.h
@@ -0,0 +1,43 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "cmath"
+
+
+namespace MuonCalib {
+
+
+
+class MeanRMS
+ {
+ public:
+ inline MeanRMS() : m_sum_x(0.), m_sum_x_sq(0.), n_ent(0) {}
+ inline void Insert(double x)
+ {
+ m_sum_x += x;
+ m_sum_x_sq += x*x;
+ n_ent++;
+ }
+ inline double GetMean() const
+ {
+ if(n_ent==0) return 0;
+ return m_sum_x/static_cast<double>(n_ent);
+ }
+ inline double GetRMS() const
+ {
+ if(n_ent==0) return 0;
+ double mean=GetMean();
+ return std::sqrt(m_sum_x_sq/static_cast<double>(n_ent) - mean*mean);
+ }
+ inline int GetN() const
+ {
+ return n_ent;
+ }
+ private:
+ double m_sum_x;
+ double m_sum_x_sq;
+ int n_ent;
+
+ };
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationAnalytic.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationAnalytic.h
new file mode 100644
index 0000000000000..83d4ac66fa308
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationAnalytic.h
@@ -0,0 +1,398 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 06.04.2006, AUTHOR: OLIVER KORTNER
+// Modified: 02.07.2006 by O. Kortner, redesign to simplified optimization.
+// 16.07.2006 by O. Kortner, optimized quality cuts.
+// 17.07.2006 by O. Kortner, compatibility with old event loop.
+// 13.12.2006 by O. Kortner and J. von Loeben, new convergence
+// criterion.
+// 27.03.2007 by O. Kortner, control histograms are written to file
+// when they are switched off again.
+// 06.08.2007 by O. Kortner, option to force r(t) to be monotonically
+// increasing
+// 29.10.2007 by O. Kortner, avoid double analysis in case of full
+// chamber track fits.
+// 30.06.2008 by O. Kortner, m_r_max settings corrected, better hit
+// selection.
+// 03.08.2008 by O. Kortner, r-t smoothing using the conventional
+// autocalibration method is implemented.
+// 28.02.2009 by O. Kortner, parabolic extrapolation added.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef RtCalibrationAnalyticH
+#define RtCalibrationAnalyticH
+
+//:::::::::::::::::::::::::::::::::
+//:: CLASS RtCalibrationAnalytic ::
+//:::::::::::::::::::::::::::::::::
+
+namespace MuonCalib {
+/// \class RtCalibrationAnalytic
+/// This class performs the analytic autocalibration whose basic ideas were
+/// developed by Mario Deile (see ATL-MUON-2004-021).
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 06.04.2006
+/**
+@class RtCalibrationAnalytic
+This class performs the analytic autocalibration whose basic ideas were
+developed by Mario Deile (see ATL-MUON-2004-021).
+
+@author Oliver.Kortner@CERN.CH
+
+*/
+
+}
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// STL //
+#include <vector>
+#include <string>
+
+// CLHEP //
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+#include "CLHEP/Matrix/SymMatrix.h"
+#include "CLHEP/Matrix/Vector.h"
+
+// ROOT //
+#include "TFile.h"
+#include "TH1F.h"
+#include "TH2F.h"
+
+// MuonCalib //
+#include "MdtCalibInterfaces/IMdtCalibration.h"
+#include "MdtCalibFitters/QuasianalyticLineReconstruction.h"
+#include "MdtCalibRt/MeanRMS.h"
+
+namespace MuonCalib {
+
+class IRtRelation;
+class RtRelationLookUp;
+class RtCalibrationOutput;
+class IMdtCalibrationOutput;
+class MuonCalibSegment;
+class BaseFunction;
+
+class RtCalibrationAnalytic : public IMdtCalibration {
+
+public:
+// Constructors //
+ RtCalibrationAnalytic(std::string name) : IMdtCalibration(name) {
+ init(0.5*CLHEP::mm, 1, 5, true, true, true, true, 100, false, false);
+ }
+ ///< 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.
+
+ RtCalibrationAnalytic(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.
+ ///< \param func_type: type of function to be used for the r-t correction;
+ ///< = 1: Legendre polynomial,
+ ///< = 2: Chebyshev polynomial,
+ ///< = 3: polygon equidistant in r.
+ ///< \param ord The order of the r-t correction polynomial is set to ord.
+ ///< \param split = true forces the algorithm to restrict segments to
+ ///< multilayers. Otherwise, segments may contain hits from different
+ ///< multilayers.
+ ///< \param 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.
+ ///< \param fix_min,fix_max =true fix r(t_min), r(t_max) (this is default).
+ ///< \param max_it maximum number of iterations.
+ ///< \param do_smoothing Smoothen the r-t relations after convergence.
+ ///> \param do_parabolic_extrapolation Use parabolic extrapolations for
+ ///< small and larged drift radii.
+
+// Desctructor //
+ ~RtCalibrationAnalytic(void);
+
+// Methods //
+// get-methods //
+ double reliability(void) 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(void) 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(void) const;
+ ///< get the number of segments which
+ ///< were passed to the algorithm
+ int numberOfSegmentsUsed(void) const;
+ ///< get the number of segments which
+ ///< are used in the autocalibration
+ int iteration(void) const;
+ ///< get the number of the current
+ ///< iteration
+ bool splitIntoMultilayers(void) const;
+ ///< returns true, if segments are
+ ///< internally restricted to single
+ ///< multilayers;
+ ///< returns false, if segments may
+ ///< contain hits from both multilayers
+ bool fullMatrix(void) 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(void) const;
+ ///< returns true, if the r-t relationship
+ ///< will be smoothened using the
+ ///< conventional autocalibration after
+ ///< convergence;
+ ///< returns false otherwise
+
+// set-method //
+ 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(void);
+ ///< the algorithm does not produce
+ ///< controll histograms (this is the
+ ///< default)
+ void forceMonotony(void);
+ ///< force r(t) to be monotonically
+ ///< increasing
+ void doNotForceMonotony(void);
+ ///< do not force r(t) to be
+ ///< monotonically increasing
+ void doSmoothing(void);
+ ///< requires that the r-t relationship
+ ///< will be smoothened using the
+ ///< conventional autocalibration after
+ ///< convergence
+ void noSmoothing(void);
+ ///< do not smoothen the r-t relationship
+ ///< after convergence
+ void doParabolicExtrapolation(void);
+ ///< requires that parabolic extrapolation
+ ///< will be used for small and large radii
+ void noParabolicExtrapolation(void);
+ ///< no parabolic extrapolation is done
+
+// methods required by the base class "IMdtCalibration" //
+ const IMdtCalibrationOutput * analyseSegments(
+ const std::vector<MuonCalibSegment*> & 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(void);
+ ///< perform the autocalibration with
+ ///< the segments acquired so far
+ bool converged(void) const;
+ ///< returns true, if the
+ ///< autocalibration has converged
+ const IMdtCalibrationOutput * getResults(void) const;
+ ///< returns the final r-t relationship
+
+private:
+// options //
+ bool m_control_histograms; // = true, if control histograms should be
+ // produces
+ bool m_split_into_ml; // = true, if segments should be restricted to the
+ // multilayers;
+ // = false, if segments over both multilayers are
+ // allowed
+ bool m_full_matrix; // = true, if the full matrix relating the errors in
+ // the r-t relationship to the residuals
+ // should be used;
+ // = false, if a diagonal matrix should be used;
+ // in this case the algorithm is equivalent
+ // to conventional method
+ bool m_fix_min; // = true: fix r(t_min)
+ bool m_fix_max; // = true: fix r(t_max)
+ int m_max_it; // maximum number of iterations
+ bool m_force_monotony; // = true if r(t) is forced to monotonically
+ // increasing, false otherwise
+
+// bookkeeping //
+ int m_nb_segments; // number of segments passed to the algorithm
+ int m_nb_segments_used; // number of segments used by the algorithm
+ int m_iteration; // current iteration
+ bool m_multilayer[2]; // m_multilayer[k] = true, if there was a segment
+ // extending to multilayer k+1
+
+// r-t quality //
+ int m_status; // m_status: 0: no covergence yet,
+ // 1: convergence, r-t is reliable,
+ // 2: convergence, r-t is unreliable
+ double m_rt_accuracy; // 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;
+ // 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; // 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 //
+ const IRtRelation * m_rt; // pointer to the input r-t relationship
+ double m_t_length; // size of the drift time interval
+ double m_t_mean; // mean value of the drift time interval
+
+// r-t output //
+ IRtRelation * m_rt_new; // r-t as determined by the autocalibration
+ RtCalibrationOutput * m_output; // class holding the results of the
+ // autocalibration
+
+// straight-segment fitting //
+ double m_r_max; // maximum value for accepted drift radii
+ QuasianalyticLineReconstruction m_tracker; // quasianalytic
+ // straight-line segment
+ // finder, used because it
+ // performs a pattern
+ // recognition
+
+// autocalibration objects //
+ bool m_do_smoothing; // = true: the r-t relationship is smoothened after
+ // convergence, no smoothing is done
+ // otherwise
+ bool m_do_parabolic_extrapolation; // = true: parabolic extrapolation is
+ // done for small and large radii
+ // = false: no parabolic extrapolation is
+ // done
+ unsigned int m_order; // order of the polynomial describing the
+ // correction to the r-t relationship
+ std::vector<CLHEP::HepVector> m_U; // vector of base function values
+ 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)
+
+// Legendre polynomial //
+ BaseFunction *m_base_function; // pointer to the base function u
+
+// control histograms //
+ TFile *m_tfile; // ROOT file
+ TH1F *m_cut_evolution; // cut evolution histogram
+ TH1F *m_nb_segment_hits; // number of hits on the segments
+ TH1F *m_CL; // confidence level distribution of the selected segments
+ TH2F *m_residuals; // residual distribution
+
+// private methods //
+ 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);
+ // 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;
+ // full_matrix = true forces the algorithm ti
+ // use the full matrix relating the errors in
+ // in r(t) to the residuals; otherwise a unit
+ // matrix is used;
+ // fix_min, fix_max=true: fix r(t_min), r(t_max);
+ // do_smoothing: smoothen the r-t relationship
+ // after convergence if do_smoothing = true;
+ // do_parabolic_extrapolation: do or not do
+ // parabolic extrapolation for small and large
+ // radii;
+ // max_it: maximum number of iterations
+ 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);
+
+ 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;
+
+ MeanRMS m_track_slope;
+ MeanRMS m_track_position; //mean and rms of track slope and position
+
+
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h
new file mode 100644
index 0000000000000..adc89c73aa378
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationCurved.h
@@ -0,0 +1,352 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 10.05.2008, AUTHOR: OLIVER KORTNER
+// Modified: 01.03.2009 by O. Kortner, parabolic extrapolation added.
+// 15.03.2009 by O. Kortner, smoothing added.
+// 18.03.2009 by O. Kortner, method performParabolicExtrapolation
+// added.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef MuonCalib_RtCalibrationCurvedH
+#define MuonCalib_RtCalibrationCurvedH
+
+//:::::::::::::::::::::::::::::::
+//:: CLASS RtCalibrationCurved ::
+//:::::::::::::::::::::::::::::::
+
+/// \class RtCalibrationCurved
+///
+/// This class performs the autocalibration of MDT chambers from the residuals
+/// of curved (parabolic) segments. The method used is a generalization of the
+/// analytic autocalibration developed by Deile, Hessey, and Staude (see
+/// ATL-MUON-2004-021).
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 10.05.2008
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// STL //
+#include <vector>
+#include <string>
+#include <list>
+
+// CLHEP //
+#include "CLHEP/Matrix/SymMatrix.h"
+#include "CLHEP/Matrix/Vector.h"
+
+// ROOT //
+#include "TFile.h"
+#include "TH1F.h"
+#include "TH2F.h"
+
+// MuonCalib //
+#include "MdtCalibInterfaces/IMdtCalibration.h"
+
+#include "GeoPrimitives/GeoPrimitives.h"
+
+namespace MuonCalib {
+
+class IMdtCalibrationOutput;
+class IRtRelation;
+class RtRelationLookUp;
+class RtCalibrationOutput;
+class CurvedPatRec;
+class MuonCalibSegment;
+class BaseFunction;
+class Legendre_polynomial;
+class CurvedLine;
+class MultilayerRtDifference;
+
+class RtCalibrationCurved : public IMdtCalibration {
+
+public:
+// Constructors //
+ RtCalibrationCurved(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.
+
+ RtCalibrationCurved(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.
+ ///< \param rt_accuracy r-t accuracy in mm. The r-t accuracy is used
+ ///< internally to distinguish between good and bad segments.
+ ///< \param 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.
+ ///< \param ord Order of the r-t correction function (15 is default).
+ ///< \param fix_min =true: r(t0) is fixed in the autocalibration procedure
+ ///< (this is default).
+ ///< \param fix_max =true: r(tmax) is fixed in the autocalibration procedure
+ ///< (false is default).
+ ///< \param max_it Maximum number of iterations (20 by default).
+ ///< \param do_parabolic_extrapolation Use parabolic extrapolations for
+ ///< small and larged drift radii.
+ ///< \param do_smoothing Smoothing of the r-t relations after convergence.
+
+// Destructor //
+ ~RtCalibrationCurved(void);
+ ///< Destructor.
+
+// Methods //
+// get-methods //
+ double reliability(void) 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(void) 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(void) const;
+ ///< get the number of segments which
+ ///< were passed to the algorithm
+ int numberOfSegmentsUsed(void) const;
+ ///< get the number of segments which
+ ///< are used in the autocalibration
+ int iteration(void) const;
+ ///< get the number of the current
+ ///< iteration
+ bool smoothing(void) const;
+ ///< returns true, if the r-t relationship
+ ///< will be smoothened using the
+ ///< conventional autocalibration after
+ ///< convergence;
+ ///< returns false otherwise
+
+// set-method //
+ 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(void);
+ ///< the algorithm does not produce
+ ///< controll histograms (this is the
+ ///< default)
+ void forceMonotony(void);
+ ///< force r(t) to be monotonically
+ ///< increasing (this is default)
+ void doNotForceMonotony(void);
+ ///< do not force r(t) to be
+ ///< monotonically increasing
+ void doParabolicExtrapolation(void);
+ ///< requires that parabolic extrapolation
+ ///< will be used for small and large radii
+ void noParabolicExtrapolation(void);
+ ///< no parabolic extrapolation is done
+ void doSmoothing(void);
+ ///< requires that the r-t relationship
+ ///< will be smoothened using the
+ ///< conventional autocalibration after
+ ///< convergence
+ void noSmoothing(void);
+ ///< do not smoothen the r-t relationship
+ ///< after convergence
+
+// methods required by the base class "IMdtCalibration" //
+ const IMdtCalibrationOutput * analyseSegments(
+ const std::vector<MuonCalibSegment*> & 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(const std::vector<MuonCalibSegment*> & seg);
+ ///< perform the autocalibration with
+ ///< the segments acquired so far
+ bool converged(void) const;
+ ///< returns true, if the
+ ///< autocalibration has converged
+ const IMdtCalibrationOutput * getResults(void) const;
+ ///< returns the final r-t relationship
+
+private:
+// options //
+ bool m_control_histograms; // = true, if control histograms should be
+ // produces
+ bool m_fix_min; // = true: fix r(t_min)
+ bool m_fix_max; // = true: fix r(t_max)
+ int m_max_it; // maximum number of iterations
+ bool m_force_monotony; // = true if r(t) is forced to monotonically
+ // increasing, false otherwise
+ bool m_do_multilayer_rt_scale; // determine multilayer rt scaling
+
+// bookkeeping //
+ int m_nb_segments; // number of segments passed to the algorithm
+ int m_nb_segments_used; // number of segments used by the algorithm
+ int m_iteration; // current iteration
+ bool m_multilayer[2]; // m_multilayer[k] = true, if there was a segment
+ // extending to multilayer k+1
+
+// r-t quality //
+ int m_status; // m_status: 0: no covergence yet,
+ // 1: convergence, r-t is reliable,
+ // 2: convergence, r-t is unreliable
+ double m_rt_accuracy; // 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;
+ // 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; // 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 //
+ const IRtRelation * m_rt; // pointer to the input r-t relationship
+ double m_t_length; // size of the drift time interval
+ double m_t_mean; // mean value of the drift time interval
+
+// r-t output //
+ IRtRelation * m_rt_new; // r-t as determined by the autocalibration
+ RtCalibrationOutput * m_output; // class holding the results of the
+ // autocalibration
+ MultilayerRtDifference * m_multilayer_rt_difference;
+// curved-segment fitting //
+ double m_r_max; // maximum value for accepted drift radii
+ 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; // = true: parabolic extrapolation is
+ // done for small and large radii
+ // = false: no parabolic extrapolation is
+ // done
+ bool m_do_smoothing; // = true: the r-t relationship is smoothened after
+ // convergence, no smoothing is done
+ // otherwise
+ unsigned int m_order; // 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 //
+ BaseFunction *m_base_function; // pointer to the base function u
+ Legendre_polynomial *m_Legendre; // pointer to the Legendre polynomial
+ // describing the curved line
+
+// control histograms //
+ TFile *m_tfile; // ROOT file
+ TH1F *m_cut_evolution; // cut evolution histogram
+ TH1F *m_nb_segment_hits; // number of hits on the segments
+ TH1F *m_pull_initial; // initial pull distribution
+ TH1F *m_pull_final; // final pull distribution after convergence
+ TH2F *m_residuals_initial; // initial residual distribution
+ TH2F *m_residuals_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
+ 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;
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationIntegration.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationIntegration.h
new file mode 100644
index 0000000000000..26cfcaaa1a5a9
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationIntegration.h
@@ -0,0 +1,127 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 02.02.2007, AUTHOR: OLIVER KORTNER
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef MuonCalib_RtCalibrationIntegrationH
+#define MuonCalib_RtCalibrationIntegrationH
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// STL //
+#include <string>
+#include <vector>
+
+// MuonCalib //
+#include "MdtCalibInterfaces/IMdtCalibration.h"
+#include "MdtCalibInterfaces/IMdtCalibrationOutput.h"
+#include "MdtCalibData/IRtRelation.h"
+#include "MdtCalibRt/RtCalibrationOutput.h"
+#include "MdtCalibInterfaces/IMdtSegmentFitter.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+
+namespace MuonCalib {
+
+//::::::::::::::::::::::::::::::::::::
+//:: CLASS RtCalibrationIntegration ::
+//::::::::::::::::::::::::::::::::::::
+
+/// \class RtCalibrationIntegration
+///
+/// This class allows the user to obtain an r-t relationship from the drift
+/// time spectrum in a given calibration region. The user can ask the class
+/// to restrict the r-t determination on hit on the segments or to use close
+/// hits too. The algorithms performs a t0 and tmax fit in order to get t(r=0)
+/// and t(r=rmax) right.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 02.02.2007
+
+class RtCalibrationIntegration : public IMdtCalibration {
+
+public:
+// Constructors //
+ RtCalibrationIntegration(std::string name) : IMdtCalibration(name) {
+ init(false, 14.6,13.0,14.0, false);
+ }
+ ///< Default constructor. Only hits on the segment are used by the
+ ///< algorithm by default. The maximum drift radius is set to 14.6 mm.
+
+ RtCalibrationIntegration(std::string name,
+ bool close_hits,
+ double r_max,
+ double lower_extrapolation_radius,
+ double upper_extrapolation_radius,
+ bool add_tmax_difference
+ ) : IMdtCalibration(name) {
+ init(close_hits, r_max, lower_extrapolation_radius, upper_extrapolation_radius, add_tmax_difference);
+ }
+ ///< Constructor. If close_hits is set to true, also close hits are
+ ///< used. The maximum drift radius is set to r_max [mm].
+
+// Methods //
+// get-methods //
+ unsigned int number_of_hits_used(void) const;
+ ///< get the number of hits used in the
+ ///< r-t determination
+
+// methods required by the base class "IMdtCalibration" //
+ const IMdtCalibrationOutput * analyseSegments(
+ const std::vector<MuonCalibSegment*> & seg);
+ ///< determine r(t)
+ bool handleSegment(MuonCalibSegment & seg);
+ ///< analyse the segment "seg"
+ void setInput(const IMdtCalibrationOutput * rt_input);
+ ///< the method is empty as no initial
+ ///< r-t relationship is required by
+ ///< the algorithm
+ bool analyse(void);
+ ///< perform the integration method
+ bool converged(void) const;
+ ///< returns true, if the integration
+ ///< method has been performed
+ const IMdtCalibrationOutput * getResults(void) const;
+ ///< returns the final r-t relationship
+
+private:
+// options //
+ bool m_close_hits; // = true, if close hits should be used,
+ // = false, if only hits on segments should be used
+ double m_lower_extrapolation_radius;///< sets the lower radius to perform the
+ double m_upper_extrapolation_radius;///< parabolic extrapolation. the parabolic fit will be
+ ///< performet between the lower extrapolation radius
+ ///< and the upper extrapolation value and be extrapolated
+ ///< till r_max
+ bool m_add_tmax_difference;
+// drift-time spectrum //
+
+ std::vector<std::pair<double, bool> > m_t_drift; // measured drift times
+// r(t) //
+ IRtRelation * m_rt; // pointer to the final r-t relationship
+ unsigned int m_nb_hits_used; // number of hits used in the algorithm
+ unsigned int m_nb_segments_used; // number of segments used
+ double m_r_max; // maximum drift radius
+ RtCalibrationOutput * m_output; // class holding the results of the
+ // autocalibration
+
+// private methods //
+ void init(bool close_hits,
+ double r_max,
+ double lower_extrapolation_radius,
+ double higher_extrapolation_radius,
+ bool add_tmax_difference );
+ // initialize the class, close hits are used,
+ // if close_hits is set to true; the maximum
+ // drift radius is set to r_max
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationOutput.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationOutput.h
new file mode 100644
index 0000000000000..5f1b089b6f16e
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtCalibrationOutput.h
@@ -0,0 +1,36 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef RTCALIBRATIONOUTPUT_H
+#define RTCALIBRATIONOUTPUT_H
+
+#include "MdtCalibInterfaces/IMdtCalibrationOutput.h"
+#include "MdtCalibData/RtFullInfo.h"
+#include "MdtCalibData/IRtRelation.h"
+
+namespace MuonCalib{
+/**
+@class RtCalibrationOutput
+Class for communication between event loop and rt calibration algorithm
+contains only a rt relation for now.
+
+*/
+
+class RtCalibrationOutput : public IMdtCalibrationOutput {
+ public:
+ RtCalibrationOutput( const IRtRelation* rt_rel, const RtFullInfo* fi ) :
+ IMdtCalibrationOutput("RtCalibrationOutput"), m_rtRelation(rt_rel), m_fullInfo(fi) {}
+
+ /** access to private attributes */
+ const IRtRelation* rt() const { return m_rtRelation; }
+ const RtFullInfo* fullInfo() const { return m_fullInfo; }
+
+ private:
+ // pointer to a IRtRelation instance
+ const IRtRelation* m_rtRelation;
+ /** additonal info for validation */
+ const RtFullInfo* m_fullInfo;
+};
+}
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtParabolicExtrapolation.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtParabolicExtrapolation.h
new file mode 100644
index 0000000000000..d367733080e89
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/MdtCalibRt/RtParabolicExtrapolation.h
@@ -0,0 +1,86 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 10.11.2008, AUTHOR: OLIVER KORTNER
+// Modified: 28.02.2009 by O. Kortner and J. von Loeben, extend extrapolation
+// features
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::::::::::::::::::::
+//:: CLASS RtParabolicExtrapolation ::
+//::::::::::::::::::::::::::::::::::::
+
+/// \class RtParabolicExtrapolation
+///
+/// This class is used to fit a parabola to the r(t) curve for in
+/// $[r_{min}, r_{max}]$ and to set r(t) according to this parabola for
+/// $r_{ext}>r>r_{max}$, if $r_{ext}>r_{max}$, and $r_{ext}<r<r_{min}$,
+/// if $r_{ext}<r_{min}$. The boundaries of the fit interval can be changed by
+/// the user.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 28.02.2009
+
+#ifndef MuonCalib_RtParabolicExtrapolationH
+#define MuonCalib_RtParabolicExtrapolationH
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+
+namespace MuonCalib {
+class IRtRelation;
+class RtRelationLookUp;
+
+class RtParabolicExtrapolation {
+
+public:
+// Constructors //
+ RtParabolicExtrapolation(void);
+ ///< Default constructor.
+
+// Methods //
+ RtRelationLookUp getRtWithParabolicExtrapolation(
+ const IRtRelation & in_rt,
+ const double & r_min=13.0,
+ const double & r_max=14.0) const;
+ ///< get an r-t relationship which is
+ ///< equivalent to the input relationship
+ ///< in_rt for r<r_min and has r(t) for
+ ///< r>r_max according to a parabola fitted
+ ///< in [r_min, r_max];
+ ///< this method is there for backward
+ ///< compatibility
+ RtRelationLookUp getRtWithParabolicExtrapolation(
+ const IRtRelation & in_rt,
+ const double & r_min,
+ const double & r_max,
+ const double & r_ext,
+ const std::vector<SamplePoint> & add_fit_points) const;
+ ///< The method fits a parabola to the r-t
+ ///< points in $[r_{min}, r_{max}]$ and the
+ ///< additional user points add_fit_points.
+ ///< The input r-t relationship in_rt is
+ ///< replaced by the parabola in
+ ///< $[r_{ext}, r_{min}]$ is
+ ///< $r_{ext}<r_{min}$ and in
+ ///< $[r_{max}, r_{ext}]$ if
+ ///< $r_{ext}>r_{max}$.
+
+private:
+ double t_from_r(const double & r, const IRtRelation & in_rt) const;
+ // get t(r) as defined in in_rt
+ double get_max_t_at_r(const double & r, const IRtRelation & in_rt) const;
+ // get the largest time of the r-t
+ // relationship in_rt corresponding to
+ // the given radius r
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/cmt/requirements b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/cmt/requirements
new file mode 100644
index 0000000000000..e2044bc0f304b
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/cmt/requirements
@@ -0,0 +1,31 @@
+package MdtCalibRt
+
+author Fabrizio Petrucci
+
+use AtlasPolicy AtlasPolicy-*
+private
+apply_tag ROOTMathLibs
+apply_tag ROOTGraphicsLibs
+apply_tag ROOTRooFitLibs
+end_private
+
+use AtlasROOT AtlasROOT-* External
+use AtlasCLHEP AtlasCLHEP-* External
+
+use MuonCalibEventBase * MuonSpectrometer/MuonCalib
+use MuonCalibMath * MuonSpectrometer/MuonCalib/MuonCalibUtils
+use MdtCalibInterfaces * MuonSpectrometer/MuonCalib/MdtCalib
+use MdtCalibData * MuonSpectrometer/MuonCalib/MdtCalib
+use MdtCalibFitters * MuonSpectrometer/MuonCalib/MdtCalib
+use GeoPrimitives * DetectorDescription
+
+include_dirs "$(MdtCalibRt_root)"
+
+library MdtCalibRt ../src/*.cxx
+apply_pattern installed_library
+
+private
+use MdtCalibT0 * MuonSpectrometer/MuonCalib/MdtCalib
+#use MuonCalibStl * MuonSpectrometer/MuonCalib/MuonCalibUtils
+#use MdtCalibUtils * MuonSpectrometer/MuonCalib/MdtCalib
+
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/doc/mainpage.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/doc/mainpage.h
new file mode 100644
index 0000000000000..013387dc5accd
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/doc/mainpage.h
@@ -0,0 +1,40 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+@mainpage MdtCalibRt Package
+@author Fabrizio.Petrucci@cern.ch
+
+@section MdtCalibRtIntro Introduction
+This package is foreseen to collect the various implementations of Mdt Rt calibration.
+
+@section MdtCalibRtOverview Class Overview
+The MdtCalibRt package contains at the moment the following classes:
+
+- MuonCalib::RtCalibrationAnalytic: This class performs the analytic autocalibration whose basic ideas were developed by Mario Deile (see ATL-MUON-2004-021)
+
+- MuonCalib::RtCalibrationIntegration: This class allows the user to obtain an r-t relationship from the drift time spectrum in a given calibration region.
+
+- MuonCalib::RtCalibrationClassic: Implementation of a rt calibration using the classical approach (iterative residuals minimization)
+- MuonCalib::ClassicSettings: Parameters needed by the Rt classic calibrations which should be provided to the RtCalibrationClassic constructor
+
+- MuonCalib::RtCalibrationOutput: Class for communication between event loop and rt calibration algorithm contains only a rt relation for now
+
+- MuonCalib::RtHistos: Set of histograms used in an iteration of RtCalibrationClassic
+
+@ref used_MdtCalibRt
+
+@ref requirements_MdtCalibRt
+*/
+
+/**
+@page used_MdtCalibRt Used Packages
+@htmlinclude used_packages.html
+*/
+
+/**
+@page requirements_MdtCalibRt Requirements
+@include requirements
+*/
+
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/AdaptiveResidualSmoothing.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/AdaptiveResidualSmoothing.cxx
new file mode 100644
index 0000000000000..5427fb0c89711
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/AdaptiveResidualSmoothing.cxx
@@ -0,0 +1,484 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 29.07.2008, AUTHOR: OLIVER KORTNER
+// 05.02.2010, JOERG v.LOEBEN
+// - tuning the method on cosmics
+// - new way to determine the binning according to satatistics
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// standard C++ //
+#include <algorithm>
+#include <iostream>
+#include <fstream>
+
+// ROOT //
+#include "TSpline.h"
+
+// MuonCalib //
+#include "MuonCalibMath/BaseFunctionFitter.h"
+#include "MuonCalibMath/PolygonBase.h"
+#include "MuonCalibMath/ConstantContentBinMaker.h"
+#include "MdtCalibRt/AdaptiveResidualSmoothing.h"
+#include "MdtCalibData/IRtRelation.h"
+#include "MuonCalibMath/DataPoint.h"
+#include "CLHEP/Matrix/Vector.h"
+
+//::::::::::::::::::::::::
+//:: NAMESPACE SETTINGS ::
+//::::::::::::::::::::::::
+
+using namespace std;
+using namespace CLHEP;
+using namespace MuonCalib;
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS ::
+//:: AdaptiveResidualSmoothing ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: CONSTRUCTOR ::
+//:::::::::::::::::
+
+AdaptiveResidualSmoothing::AdaptiveResidualSmoothing(void) {
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::
+//:: METHOD clear ::
+//::::::::::::::::::
+
+void AdaptiveResidualSmoothing::clear(void) {
+
+ m_residual_point.clear();
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD addResidual ::
+//::::::::::::::::::::::::
+
+void AdaptiveResidualSmoothing::addResidual(const double & radius,
+ const double & residual) {
+
+// make the data point //
+ CLHEP::HepVector point(2);
+ point[0] = radius;
+ point[1] = residual;
+ DataPoint residual_point(point, 0);
+
+// store the data point //
+ m_residual_point.push_back(residual_point);
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::
+//:: METHOD addResidualsFromSegment ::
+//::::::::::::::::::::::::::::::::::::
+
+bool AdaptiveResidualSmoothing::addResidualsFromSegment(
+ MuonCalibSegment & seg, bool curved, double road_width) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ CLHEP::HepVector point(2); // auxiliary point
+ IMdtPatRecFitter *fitter(0);
+
+////////////////////////////////////////
+// RESIDUAL DETERMINATION AND STORAGE //
+////////////////////////////////////////
+
+// perform fit //
+ if (curved) {
+
+ // set road width /
+ m_cfitter.setRoadWidth(road_width);
+
+ // set fitter pointer //
+ fitter = &m_cfitter;
+
+ } else {
+
+ // set road width /
+ m_sfitter.setRoadWidth(road_width);
+
+ // set fitter pointer //
+ fitter = &m_sfitter;
+
+ }
+
+ fitter->SetRefineSegmentFlag(false);
+
+ // prepare hit selection to exclude the first hit from the fit //
+// vector<unsigned int> hit_selection(seg.mdtHitsOnTrack(), 0);
+// hit_selection[0] = 1;
+
+ // perform the track fit; return in case of failure //
+ if (!fitter->fit(seg)) {
+// cout << seg.mdtHitsOnTrack() << endl;
+ return false;
+ }
+
+ // reject bad fits //
+ if (curved) {
+ if (m_cfitter.chi2PerDegreesOfFreedom()>50) {
+ return false;
+ }
+ } else {
+ if (m_sfitter.chi2PerDegreesOfFreedom()>50) {
+ return false;
+ }
+ }
+
+ // store the residuals //
+// static ofstream resfile("res.txt");
+ for (unsigned int k=0; k<fitter->trackHits().size(); k++) {
+ point[0] = fitter->trackHits()[k]->driftRadius();
+ point[1] = fitter->trackHits()[k]->driftRadius()-
+ fabs(fitter->trackHits()[k]->signedDistanceToTrack());
+// resfile << point[0] << "\t" << point[1] << endl;
+ DataPoint residual_point(point, 0);
+ m_residual_point.push_back(residual_point);
+ }
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::::::
+//:: METHOD performSmoothing(., ., ., .) ::
+//:::::::::::::::::::::::::::::::::::::::::
+
+RtRelationLookUp AdaptiveResidualSmoothing::performSmoothing(
+ const IRtRelation & rt_rel,
+ unsigned int nb_entries_per_bin,
+ bool fix_t_min, bool fix_t_max) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ ConstantContentBinMaker bin_maker(m_residual_point, 0.001);
+ vector<unsigned int> ref_coord(1);
+
+////////////
+// CHECKS //
+////////////
+
+ if (m_residual_point.size()==0) {
+ cerr << endl
+ << "Class AdaptiveResidualSmoothing, method performSmoothing: "
+ << "ERROR!\n"
+ << "No residuals stored.\n";
+ exit(1);
+ }
+
+ if (m_residual_point.size()<nb_entries_per_bin) {
+ cerr << endl
+ << "Class AdaptiveResidualSmoothing, method performSmoothing: "
+ << "ERROR!\n"
+ << "Not enough residuals stored.\n";
+ exit(1);
+ }
+
+/////////////////////
+// PERFORM BINNING //
+/////////////////////
+
+ ref_coord[0] = 0;
+ bin_maker.binDataPoints(nb_entries_per_bin, ref_coord);
+
+//////////////////////////////
+// DETERMINE r-t COORECTION //
+//////////////////////////////
+
+// get a polygon through the bin centres //
+// cout << "bin_maker.getBins().size()="
+// << bin_maker.getBins().size() << endl;
+ vector<double> rad(bin_maker.getBins().size());
+ vector<SamplePoint> corr(bin_maker.getBins().size());
+ for (unsigned int k=0; k<bin_maker.getBins().size(); k++) {
+ rad[k] = (bin_maker.getBins()[k])->centreOfGravity()[0];
+// cout << "rad[" << k << "] = " << rad[k] << endl;
+ corr[k].set_x1(rad[k]);
+ corr[k].set_x2((bin_maker.getBins()[k])->centreOfGravity()[1]);
+ corr[k].set_error(1.0);
+ }
+ sort(rad.begin(), rad.end());
+
+ vector<double> radi;
+ radi.push_back(rad[0]);
+ unsigned int counter(0);
+ for (unsigned int k=1; k<rad.size(); k++) {
+ if (rad[counter]<rad[k]) {
+ radi.push_back(rad[k]);
+ counter++;
+ }
+ }
+
+ PolygonBase polygon(radi);
+ BaseFunctionFitter fitter(bin_maker.getBins().size());
+ fitter.fit_parameters(corr, 1, corr.size(), &polygon);
+
+// create an improved r-t relationship //
+// vector<double> rt_point;
+// unsigned int nb_rt_points(60);
+// double step_size((rt_rel.radius(rt_rel.tUpper())-
+// rt_rel.radius(rt_rel.tLower()))/
+// static_cast<double>(nb_rt_points));
+// for (double r=rt_rel.radius(rt_rel.tLower());
+// r<=rt_rel.radius(rt_rel.tUpper()); r=r+step_size) {
+// double t(t_from_r(rt_rel, r));
+// rt_point.push_back(t);
+// double delta_r(0.0);
+// for (unsigned int l=0; l<radi.size(); l++) {
+// delta_r = delta_r+fitter.coefficients()[l]*polygon.value(l,
+// rt_rel.radius(t));
+// }
+// rt_point.push_back(rt_rel.radius(t)-delta_r);
+//
+// }
+// RtSpline improved_rt(rt_point);
+ vector<double> rt_params;
+ rt_params.push_back(rt_rel.tLower());
+ rt_params.push_back(0.01*(rt_rel.tUpper()-rt_rel.tLower()));
+// ofstream outfile("out2.txt");
+ for (double t=rt_rel.tLower(); t<=rt_rel.tUpper(); t=t+rt_params[1]) {
+ double delta_r(0.0);
+ for (unsigned int l=0; l<radi.size(); l++) {
+ delta_r = delta_r+fitter.coefficients()[l]*polygon.value(l,
+ rt_rel.radius(t));
+ }
+ if (fix_t_min && (rt_rel.radius(t)<0.5)) {
+ delta_r = 0.0;
+ }
+ if (fix_t_max && (rt_rel.radius(t)>14.1)) {
+ delta_r = 0.0;
+ }
+ rt_params.push_back(rt_rel.radius(t)-delta_r);
+ }
+
+ RtRelationLookUp improved_rt(rt_params);
+// for (double t=rt_rel.tLower(); t<=rt_rel.tUpper(); t=t+rt_params[1]) {
+// outfile << rt_rel.radius(t) << "\t" << improved_rt.radius(t)
+// << "\t" << endl;
+// }
+
+ return improved_rt;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::::
+//:: METHOD performSmoothing(., ., .) ::
+//::::::::::::::::::::::::::::::::::::::
+
+RtRelationLookUp AdaptiveResidualSmoothing::performSmoothing(
+ const IRtRelation & rt_rel,
+ const bool & fix_t_min,
+ const bool & fix_t_max) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ unsigned int nb_bins; // number of bins in r
+ unsigned int min_nb_entries_per_bin(20); // minimum number of
+ // entries per bin
+ double step_size; // real step size [mm]
+ double r_max(16.0); // maximum drift radius
+ double aux_rad; // auxiliary radius
+ double aux_corr; // auxiliary correction value
+ unsigned int bin_index;
+
+//////////////////////////////////
+// DETERMINE THE NUMBER OF BINS //
+//////////////////////////////////
+
+
+ // take the sqrt from the number of residuals to determine the binsize.
+ // for low stat, a larger statistical error on the bins is allowed
+ // for example:
+ // 1000 segments -> 7.6% and 35 bins
+ // 5000 segments -> 5.1% and 64 bins
+ // 10000 segments-> 4.1% and 91 bins
+ nb_bins =static_cast<unsigned int>( sqrt(m_residual_point.size()/6.) );
+
+ //calculate min_nb_per_bin
+ min_nb_entries_per_bin = m_residual_point.size() / nb_bins;
+ //cout << "min number of entries/bin: " << min_nb_entries_per_bin << endl;
+
+//////////////////////////////////////
+// RETURN IF THERE ARE NO RESIDUALS //
+//////////////////////////////////////
+
+ if (m_residual_point.size()==0) {
+ cerr << endl
+ << "Class AdaptiveResidualSmoothing, "
+ << "performSmoothing(., ., .): WARNING!\n"
+ << "No residuals are stored. no correction applied to r(t)!\n"
+ << endl;
+ }
+
+ if (nb_bins==0) {
+ nb_bins = 1;
+ }
+ step_size = r_max/static_cast<double>(nb_bins);
+
+//////////////////////////////////////////////////////////////
+// DETERMINE THE AVERAGE RESIDUAL VALUES IN THE RADIAL BINS //
+//////////////////////////////////////////////////////////////
+
+// vector for the correction function //
+ vector<SamplePoint> corr(nb_bins);
+ vector<double> radii(nb_bins);
+
+// sort the residuals point in the residual value for a simple outlyer //
+// rejection performed later //
+ for (unsigned int k=0; k<m_residual_point.size(); k++) {
+ m_residual_point[k].setReferenceComponent(1);
+ }
+ sort(m_residual_point.begin(), m_residual_point.end());
+
+// auxiliary data point arrays //
+ vector< vector<const DataPoint *> > sample_points_per_r_bin(nb_bins);
+
+// group data points in r bins //
+ for (unsigned int k=0; k<m_residual_point.size(); k++) {
+ if (fabs(m_residual_point[k].dataVector()[0])>=r_max) {
+ continue;
+ }
+ bin_index = static_cast<unsigned int>(fabs(
+ m_residual_point[k].dataVector()[0])/step_size);
+ sample_points_per_r_bin[bin_index].push_back(&(m_residual_point[k]));
+ }
+
+// loop over the r bins and determine the r-t corrections //
+// static ofstream outfile("tmp.txt");
+// outfile << endl;
+ for (unsigned int bin=0; bin<nb_bins; bin++) {
+
+ aux_rad = 0.0;
+ aux_corr = 0.0;
+
+ // remove the first and last 1% of the residuals as simple outlyer removal //
+ unsigned int k_min(static_cast<unsigned int>(
+ 0.01*sample_points_per_r_bin[bin].size()));
+ unsigned int k_max(static_cast<unsigned int>(
+ 0.99*sample_points_per_r_bin[bin].size()));
+// outfile << bin << "\t" << (0.5+bin)*step_size << "\t"
+// << k_min << "\t" << k_max << endl;
+
+ // loop over residuals in a bin
+ for (unsigned int k=k_min; k<k_max; k++) {
+ aux_rad = aux_rad+
+ sample_points_per_r_bin[bin][k]->dataVector()[0];
+ aux_corr = aux_corr+
+ sample_points_per_r_bin[bin][k]->dataVector()[1];
+ }
+ if (k_max-k_min<0.25*min_nb_entries_per_bin) {
+ radii[bin] = (0.5+bin)*step_size;
+ corr[bin] = SamplePoint(radii[bin], 0.0, 1.0);
+ } else {
+ radii[bin] = aux_rad/static_cast<double>(k_max-k_min);
+ corr[bin] = SamplePoint(radii[bin],
+ aux_corr/static_cast<double>(k_max-k_min), 1.0);
+ }
+
+ }
+
+// correction polygon //
+ PolygonBase polygon(radii);
+ BaseFunctionFitter fitter(nb_bins);
+ fitter.fit_parameters(corr, 1, corr.size(), &polygon);
+
+// create output r-t relationship //
+ vector<double> rt_params;
+ rt_params.push_back(rt_rel.tLower());
+ rt_params.push_back(0.01*(rt_rel.tUpper()-rt_rel.tLower()));
+// ofstream outfile("out2.txt");
+ for (double t=rt_rel.tLower(); t<=rt_rel.tUpper(); t=t+rt_params[1]) {
+ double delta_r(0.0);
+ for (unsigned int l=0; l<radii.size(); l++) {
+ delta_r = delta_r+fitter.coefficients()[l]*polygon.value(l,
+ rt_rel.radius(t));
+ }
+ if (rt_rel.radius(t)>r_max) {
+ delta_r = 0.0;
+ }
+ if (fix_t_min && (rt_rel.radius(t)<0.5)) {
+ delta_r = 0.0;
+ }
+ if (fix_t_max && (rt_rel.radius(t)>14.1)) {
+ delta_r = 0.0;
+ }
+ rt_params.push_back(rt_rel.radius(t)-delta_r);
+ }
+
+ RtRelationLookUp improved_rt(rt_params);
+// for (double t=rt_rel.tLower(); t<=rt_rel.tUpper(); t=t+rt_params[1]) {
+// outfile << rt_rel.radius(t) << "\t" << improved_rt.radius(t)
+// << "\t" << endl;
+// }
+
+ return improved_rt;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD t_from_r ::
+//:::::::::::::::::::::
+
+double AdaptiveResidualSmoothing::t_from_r(const IRtRelation & rt_rel,
+ const double & r) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ double precision(0.001); // spatial precision of the inversion
+ double t_max(rt_rel.tUpper()); // upper time search limit
+ double t_min(rt_rel.tLower()); // lower time search limit
+
+/////////////////////////////////////////////
+// SEARCH FOR THE CORRESPONDING DRIFT TIME //
+/////////////////////////////////////////////
+
+ while (t_max-t_min>0.1 &&
+ fabs(rt_rel.radius(0.5*(t_min+t_max))-r)>precision) {
+
+ if (rt_rel.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);
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.cxx
new file mode 100644
index 0000000000000..9dac19300b028
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.cxx
@@ -0,0 +1,220 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MultilayerRtDifference.h"
+
+//MuonCalibEventBase
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+
+//MdtCalibFitters
+#include "MdtCalibFitters/CurvedPatRec.h"
+#include "MdtCalibFitters/QuasianalyticLineReconstruction.h"
+
+//MdtCalibData
+#include "MdtCalibData/RtScaleFunction.h"
+#include "MdtCalibData/IRtRelation.h"
+
+
+//root
+#include "TH2F.h"
+#include "TProfile.h"
+#include "TF2.h"
+#include "TF1.h"
+#include "TDirectory.h"
+
+//std
+#include "cmath"
+#include "sstream"
+#include "iostream"
+
+namespace MuonCalib {
+
+Double_t MultilayerRtDifference_ScaleFunction(Double_t *x, Double_t *par)
+ {
+ return par[0] * RtScalePolynomial(x[0]);
+ }
+
+
+class MultilayerRtDifference_Histograms
+ {
+ public:
+ inline MultilayerRtDifference_Histograms(TDirectory *control_histogram_dir);
+ inline ~MultilayerRtDifference_Histograms();
+ inline TH2F * GetResHist(const int &ml);
+ inline TProfile * GetProfileDiff(const int & min_number_of_entries);
+ private:
+ std::list<TH2F *> m_residuals[2];
+ TH2F * m_current_residuals[2];
+ std::list<TProfile *> m_res_prov[2];
+ std::list<TProfile *> m_res_prov_diff;
+ TDirectory *m_control_histogram_dir;
+ };
+
+
+inline MultilayerRtDifference_Histograms::MultilayerRtDifference_Histograms(TDirectory *control_histogram_dir): m_control_histogram_dir(control_histogram_dir)
+ {
+ for(unsigned int i=0; i<2; i++)
+ {
+ m_current_residuals[i]=NULL;
+ }
+ }
+
+
+template<typename T> inline void clearall(T & container)
+ {
+ typename T::iterator it=container.begin();
+ for(; it!=container.end(); it++)
+ {
+ delete *it;
+ *it=NULL;
+ }
+ }
+
+inline MultilayerRtDifference_Histograms::~MultilayerRtDifference_Histograms()
+ {
+ for(unsigned int i=0; i<2; i++)
+ {
+ clearall(m_residuals[i]);
+ clearall(m_res_prov[i]);
+ }
+ clearall(m_res_prov_diff);
+ }
+
+
+inline TH2F * MultilayerRtDifference_Histograms::GetResHist(const int &ml)
+ {
+ if(!m_current_residuals[ml])
+ {
+ for(int i=0; i<2; i++)
+ {
+ TDirectory *prevdir=gDirectory;
+ if (m_control_histogram_dir)
+ m_control_histogram_dir->cd();
+ std::ostringstream nm;
+ nm << "temporal_residual_ml" << i << "_iteration" << m_residuals[i].size();
+ m_current_residuals[i]=new TH2F(nm.str().c_str(), "", 15, 0., 15., 200, -100, 100);
+ if (m_control_histogram_dir)
+ prevdir->cd();
+ else
+ m_current_residuals[i]->SetDirectory(0);
+ m_residuals[i].push_back(m_current_residuals[i]);
+ }
+ }
+ return m_current_residuals[ml];
+ }
+
+inline TProfile * MultilayerRtDifference_Histograms::GetProfileDiff(const int & min_number_of_entries)
+ {
+ TDirectory *prev=gDirectory;
+ if(m_control_histogram_dir)
+ {
+ m_control_histogram_dir->cd();
+ }
+ bool ok(true);
+ TProfile *prof[2];
+ for(unsigned int i=0; i<2; i++)
+ {
+ if(m_current_residuals[i] && m_current_residuals[i]->GetEntries()>=min_number_of_entries)
+ {
+ prof[i]=m_current_residuals[i]->ProfileX();
+ if(!m_control_histogram_dir)
+ {
+ prof[i]->SetDirectory(0);
+ }
+ }
+ else
+ {
+ prof[i]=NULL;
+ ok=false;
+ }
+ m_res_prov[i].push_back(prof[i]);
+ m_current_residuals[i]=NULL;
+ }
+ if(!ok)
+ {
+ prev->cd();
+ m_res_prov_diff.push_back(NULL);
+ return NULL;
+ }
+ std::ostringstream nm;
+ nm<<"res_prov_diff_step"<<m_res_prov_diff.size();
+ TProfile *prov_diff=new TProfile(nm.str().c_str(), "", 15, 0., 15.);
+ prov_diff->Add(prof[0], prof[1], 1., -1);
+ m_res_prov_diff.push_back(prov_diff);
+ if(!m_control_histogram_dir)
+ {
+ prov_diff->SetDirectory(0);
+ }
+ return prov_diff;
+ }
+
+MultilayerRtDifference :: MultilayerRtDifference(int min_hits, TDirectory *control_histogram_dir): m_histograms(new MultilayerRtDifference_Histograms(control_histogram_dir)), m_min_number_of_hits(min_hits)
+ {
+ m_polfun = new TF1("polfun", MultilayerRtDifference_ScaleFunction, 4.0, 15.0, 1);
+ }
+
+
+MultilayerRtDifference :: ~MultilayerRtDifference()
+ {
+ delete m_histograms;
+ delete m_polfun;
+ }
+
+void MultilayerRtDifference :: Fill(const MdtCalibHitBase &hit, const IRtRelation &rt_relation)
+ {
+ int ml=hit.identify().mdtMultilayer() -1;
+ double r_track=std::fabs(hit.signedDistanceToTrack());
+ double res=std::fabs(hit.driftRadius()) - r_track;
+ double v=rt_relation.driftvelocity(hit.driftTime());
+ m_histograms->GetResHist(ml)->Fill(r_track, res/v);
+ }
+
+bool MultilayerRtDifference :: DoFit(IRtRelation * rt_relation, const std::vector<MuonCalibSegment *> *seg)
+ {
+ TProfile *prov_diff=m_histograms->GetProfileDiff(m_min_number_of_hits);
+ if(!prov_diff)
+ {
+ std::cout<<"MultilayerRtDifference :: DoFit: Not enough hits!";
+ return false;
+ }
+ if(prov_diff->Fit("polfun", "Q", "", 4., 15.)!=0)
+ {
+ std::cerr<<"MultilayerRtDifference: Fit of polinomial failed! Not updating scale!"<<std::endl;
+ return false;
+ }
+ if(std::fabs(m_polfun->GetParameter(0)) < m_polfun->GetParError(0))
+ {
+ std::cout<<"MultilayerRtDifference: No Scale update needed!"<<std::endl;
+ std::cout<<"Scale correction: "<<m_polfun->GetParameter(0)<<" +- "<< m_polfun->GetParError(0) <<std::endl;
+ return true;
+ }
+ if(!rt_relation) return true;
+ float scale=m_polfun->GetParameter(0);
+// std::cout<<"Scale correction: "<<m_polfun->GetParameter(0)<<" +- " <<m_polfun->GetParError(0)*std::sqrt(m_polfun->GetChisquare()/m_polfun->GetNDF())<<std::endl;
+ if(std::fabs(scale)>2) scale*=4;
+ if(rt_relation->HasTmaxDiff())
+ {
+ scale+=rt_relation->GetTmaxDiff();
+ }
+ rt_relation->SetTmaxDiff(scale);
+ if(!seg)
+ return true;
+//update input segments
+ for(std::vector<MuonCalibSegment *>::const_iterator it=seg->begin(); it!=seg->end(); it++)
+ {
+ MuonCalibSegment *seg=*it;
+ for (MuonCalibSegment::MdtHitVec::iterator h_it=seg->mdtHOTBegin (); h_it!=seg->mdtHOTEnd(); h_it++)
+ {
+ MdtCalibHitBase *hit=*h_it;
+ float old_corr=hit->TemperatureTime();
+ float corr=RtScaleFunction(hit->driftTime(), hit->identify().mdtMultilayer()==2, *rt_relation);
+ hit->setTemperatureTime(corr);
+ hit->setDriftTime(hit->driftTime() - corr + old_corr);
+ hit->setDriftRadius( rt_relation->radius(hit->driftTime()), hit->sigmaDriftRadius());
+ }
+ }
+ return true;
+ }
+
+}//namespace MuonCalib
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.h
new file mode 100644
index 0000000000000..955b568819bb0
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/MultilayerRtDifference.h
@@ -0,0 +1,47 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MdtCalibRt_MultilayerRtDifference_h
+#define MdtCalibRt_MultilayerRtDifference_h
+
+class TF1;
+class TDirectory;
+
+#include <cstddef>
+#include "vector"
+#include "list"
+
+namespace MuonCalib {
+
+class MuonCalibSegment;
+class IRtRelation;
+class MdtCalibHitBase;
+class MultilayerRtDifference_Histograms;
+
+class MultilayerRtDifference {
+
+ public:
+ MultilayerRtDifference(int min_hits, TDirectory *control_histogram_dir=NULL);
+ virtual ~MultilayerRtDifference();
+
+ void Fill(const MdtCalibHitBase &hit, const IRtRelation &rt_relation);
+
+ bool DoFit(IRtRelation * rt_relation=NULL, const std::vector<MuonCalibSegment *> *seg=NULL);
+
+ inline const TF1 * GetFunction() const
+ {
+ return m_polfun;
+ }
+
+ private:
+
+ TF1 *m_polfun;
+ MultilayerRtDifference_Histograms *m_histograms;
+ int m_min_number_of_hits;
+};
+
+
+}//namespace MuonCalib
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationAnalytic.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationAnalytic.cxx
new file mode 100644
index 0000000000000..e161c97058182
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationAnalytic.cxx
@@ -0,0 +1,1641 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 06.04.2006, AUTHOR: OLIVER KORTNER
+// Modified: 02.07.2006 by O. Kortner, redesign to simplified optimization.
+// 16.07.2006 by O. Kortner, optimized quality cuts.
+// 17.07.2006 by O. Kortner, compatibility with old event loop.
+// 13.12.2006 by O. Kortner and J. von Loeben, new convergence
+// criterion.
+// 27.03.2007 by O. Kortner, control histograms are written to file
+// when they are switched off again.
+// 06.08.2007 by O. Kortner, option to force r(t) to be monotonically
+// increasing.
+// 29.10.2007 by O. Kortner, avoid double analysis in case of full
+// chamber track fits.
+// 30.06.2008 by O. Kortner, m_r_max settings corrected, better hit
+// selection.
+// 03.08.2008 by O. Kortner, r-t smoothing using the conventional
+// autocalibration method is implemented.
+// 28.02.2009 by O. Kortner, parabolic extrapolation added.
+// 08.02.2010 by J.v.Loeben, parabolic extrapolation to end-point
+// is now done after convergence/smoothing
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "cmath"
+
+#include "MdtCalibRt/RtCalibrationAnalytic.h"
+#include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
+#include "MdtCalibRt/RtCalibrationOutput.h"
+#include "MdtCalibData/RtChebyshev.h"
+#include "MdtCalibData/RtRelationLookUp.h"
+#include "MuonCalibMath/BaseFunctionFitter.h"
+#include "MuonCalibMath/ChebyshevPolynomial.h"
+#include "MuonCalibMath/LegendrePolynomial.h"
+#include "MuonCalibMath/PolygonBase.h"
+#include "MdtCalibRt/AdaptiveResidualSmoothing.h"
+#include "MdtCalibRt/RtParabolicExtrapolation.h"
+#include "MdtCalibData/RtFromPoints.h"
+
+#include "MdtCalibData/IRtRelation.h"
+#include "MdtCalibData/RtRelationLookUp.h"
+#include "MdtCalibRt/RtCalibrationOutput.h"
+#include "MdtCalibInterfaces/IMdtCalibrationOutput.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "MuonCalibMath/BaseFunction.h"
+
+
+//root
+#include "TGraphErrors.h"
+
+//std
+#include "sstream"
+#include "cmath"
+
+//:::::::::::::::::::::::
+//:: NAMESPACE SETTING ::
+//:::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS RtCalibrationAnalytic ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//*****************************************************************************
+
+
+
+RtCalibrationAnalytic :: RtCalibrationAnalytic(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,
+ bool do_parabolic_extrapolation) : IMdtCalibration(name), m_rt(NULL) {
+ init(rt_accuracy, func_type, ord, split, full_matrix,
+ fix_min, fix_max, max_it, do_smoothing, do_parabolic_extrapolation);
+ }
+
+
+
+
+
+
+
+RtCalibrationAnalytic :: ~RtCalibrationAnalytic(void)
+ {
+ if (m_base_function!=0) {
+ delete m_base_function;
+ }
+ if (m_tfile!=0) {
+ m_tfile->Write();
+ }
+ }
+
+
+//:::::::::::::::::
+//:: METHOD 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) {
+
+/////////////////////////////
+// RESET PRIVATE VARIABLES //
+/////////////////////////////
+
+ m_r_max = 15.0*CLHEP::mm;
+ m_control_histograms = false;
+ m_tfile = 0;
+ m_cut_evolution = 0;
+ m_nb_segment_hits = 0;
+ m_CL = 0;
+ m_residuals = 0;
+ 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 = fabs(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 = abs(max_it);
+
+ if (m_order==0) {
+ cerr << "\n"
+ << "Class RtCalibrationAnalytic, method init: ERROR!"
+ << "\nOrder of the correction polynomial must be >0!\n";
+ exit(1);
+ }
+
+ 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_rt_new = 0;
+ m_output = 0;
+
+ m_U = 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 = 0;
+ cerr << "\n"
+ << "Class RtCalibrationAnalytic, method init: "
+ << "ERROR!\n"
+ << "Illegal correction function type!\n";
+ exit(1);
+ }
+ switch(func_type) {
+ case 1:
+ m_base_function = new LegendrePolynomial;
+ break;
+ case 2:
+ m_base_function = new ChebyshevPolynomial;
+ break;
+ case 3:
+ if (m_order<2) {
+ cerr << "\n"
+ << "Class RtCalibrationAnalytic, "
+ << "method init: ERROR!\n"
+ << "Order must be >2 for polygons! "
+ << "It is set to " << m_order
+ << "by the user.\n";
+ exit(1);
+ }
+ 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 = new 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;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD t_from_r ::
+//:::::::::::::::::::::
+
+double RtCalibrationAnalytic::t_from_r(const double & r) {
+
+///////////////
+// 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 &&
+ fabs(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);
+
+}
+
+//*****************************************************************************
+
+////////////////////////////
+// METHOD display_segment //
+////////////////////////////
+
+void RtCalibrationAnalytic::display_segment(MuonCalibSegment * segment,
+ ofstream & outfile) {
+
+///////////////
+// 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 << "NULL " << 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 << "NULL " << 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.m_x2()*(z_min-30.0)+aux_track.b_x2()
+ << " " << z_min-30.0
+ << " " << aux_track.m_x2()*(z_max+30.0)+aux_track.b_x2()
+ << " " << z_max+30.0 << "\n";
+
+// add a wait statement //
+ outfile << "WAIT\n";
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD reliability ::
+//::::::::::::::::::::::::
+
+double RtCalibrationAnalytic::reliability(void) const {
+
+ return m_status;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::
+//:: METHOD estimatedRtAccuracy ::
+//::::::::::::::::::::::::::::::::
+
+double RtCalibrationAnalytic::estimatedRtAccuracy(void) const {
+
+ return m_rt_accuracy;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::
+//:: METHOD numberOfSegments ::
+//:::::::::::::::::::::::::::::
+
+int RtCalibrationAnalytic::numberOfSegments(void) const {
+
+ return m_nb_segments;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::
+//:: METHOD numberOfSegmentsUsed ::
+//:::::::::::::::::::::::::::::::::
+
+int RtCalibrationAnalytic::numberOfSegmentsUsed(void) const {
+
+ return m_nb_segments_used;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD iteration ::
+//::::::::::::::::::::::
+
+int RtCalibrationAnalytic::iteration(void) const {
+
+ return m_iteration;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::
+//:: METHOD splitIntoMultilayers ::
+//:::::::::::::::::::::::::::::::::
+
+bool RtCalibrationAnalytic::splitIntoMultilayers(void) const {
+
+ return m_split_into_ml;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD fullMatrix ::
+//:::::::::::::::::::::::
+
+bool RtCalibrationAnalytic::fullMatrix(void) const {
+
+ return m_full_matrix;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD smoothing ::
+//:::::::::::::::::::::::
+
+bool RtCalibrationAnalytic::smoothing(void) const {
+
+ return m_do_smoothing;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::::::
+//:: METHOD switch_on_control_histograms ::
+//:::::::::::::::::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::switch_on_control_histograms(
+ const string & file_name) {
+
+/////////////////////////////////////////////
+// CREATE THE ROOT FILE AND THE HISTOGRAMS //
+/////////////////////////////////////////////
+
+ m_control_histograms = true;
+
+ m_tfile = new TFile(file_name.c_str(), "RECREATE");
+
+ m_cut_evolution = new TH1F("m_cut_evolution",
+ "CUT EVOLUTION",
+ 11, -0.5, 10.5);
+
+ m_nb_segment_hits = new TH1F("m_nb_segment_hits",
+ "NUMBER OF HITS ON THE REFITTED SEGMENTS",
+ 11, -0.5, 10.5);
+
+ m_CL = new TH1F("m_CL",
+ "CONFIDENCE LEVELS OF THE REFITTED SEGMENTS",
+ 100, 0.0, 1.0);
+
+ m_residuals = new TH2F("m_residuals",
+ "RESIDUALS OF THE REFITTED SEGMENTS",
+ 100, -0.5, 15.0, 300, -1.5, 1.5);
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD switch_off_control_histograms ::
+//::::::::::::::::::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::switch_off_control_histograms(void) {
+
+ m_control_histograms = false;
+ if (m_tfile!=0) {
+ m_tfile->Write();
+ }
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD forceMonotony ::
+//::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::forceMonotony(void) {
+
+ m_force_monotony = true;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::
+//:: METHOD doNotForceMonotony ::
+//:::::::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::doNotForceMonotony(void) {
+
+ m_force_monotony = false;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD doSmoothing ::
+//::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::doSmoothing(void) {
+
+ m_do_smoothing = true;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD noSmoothing ::
+//::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::noSmoothing(void) {
+
+ m_do_smoothing = false;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::
+//:: METHOD doParabolicExtrapolation ::
+//:::::::::::::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::doParabolicExtrapolation(void) {
+
+ m_do_parabolic_extrapolation = true;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::
+//:: METHOD noParabolicExtrapolation ::
+//:::::::::::::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::noParabolicExtrapolation(void) {
+
+ m_do_parabolic_extrapolation = false;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::
+//:: METHOD setEstimateRtAccuracy ::
+//::::::::::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::setEstimateRtAccuracy(const double & acc) {
+
+ m_rt_accuracy = fabs(acc);
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::
+//:: METHOD splitIntoMultilayers ::
+//:::::::::::::::::::::::::::::::::
+
+void RtCalibrationAnalytic::splitIntoMultilayers(const bool & yes_or_no) {
+
+ m_split_into_ml = yes_or_no;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD fullMatrix ::
+//:::::::::::::::::::::::
+
+void RtCalibrationAnalytic::fullMatrix(const bool & yes_or_no) {
+
+ m_full_matrix = yes_or_no;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::
+//:: METHOD analyseSegments ::
+//::::::::::::::::::::::::::::
+
+const IMdtCalibrationOutput * RtCalibrationAnalytic::analyseSegments(
+ const std::vector<MuonCalibSegment*> & seg) {
+
+ const IRtRelation *tmp_rt(0);
+ const IRtRelation *conv_rt(0);
+
+//////////////////////////
+// AUTOCALIBRATION LOOP //
+//////////////////////////
+
+ while (!converged()) {
+
+ for (unsigned int k=0; k<seg.size(); k++) {
+ handleSegment(*seg[k]);
+ }
+ if(!analyse())
+ {
+ for(unsigned int i=0; i<seg.size(); i++)
+ {
+ std::cout<<i<<" "<<seg[i]->direction()<<" "<<seg[i]->position()<<std::endl;
+ }
+ if(conv_rt) delete conv_rt;
+ return NULL;
+ }
+
+ const RtCalibrationOutput* rtOut = m_output;
+
+ if (!converged()) {
+ setInput(rtOut);
+ }
+ tmp_rt = rtOut->rt();
+
+ 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));
+ }
+ if (conv_rt) delete conv_rt;
+ conv_rt = new RtRelationLookUp(params);
+
+ }
+ if (conv_rt) delete conv_rt;
+// parabolic extrapolations for small radii //
+ if (m_do_parabolic_extrapolation) {
+ RtRelationLookUp *tmprt =
+ performParabolicExtrapolation(true, true, *tmp_rt);
+ delete m_output;
+ m_output = new RtCalibrationOutput(tmprt,
+ new RtFullInfo("RtCalibrationAnalyticExt", m_iteration,
+ m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+ delete tmp_rt;
+ 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);
+ unsigned int counter2(0);
+
+ // overwrite drift radii and calculate the average resolution //
+ for (unsigned int k=0; k<seg.size(); k++) {
+ if (seg[k]->mdtHitsOnTrack()<3) {
+ continue;
+ }
+ counter2++;
+ double avres(0.0);
+ for (unsigned int h=0; h<seg[k]->mdtHitsOnTrack(); h++) {
+ seg[k]->mdtHOT()[h]->setDriftRadius(tmp_rt->radius(
+ seg[k]->mdtHOT()[h]->driftTime()),
+ seg[k]->mdtHOT()[h]->sigmaDriftRadius());
+ if (seg[k]->mdtHOT()[h]->sigmaDriftRadius()<0.5*m_r_max) {
+ avres = avres+seg[k]->mdtHOT()[h]->sigma2DriftRadius();
+ } else {
+ avres = avres+0.1;
+ }
+ }
+ avres = avres/static_cast<double>(seg[k]->mdtHitsOnTrack());
+ avres = sqrt(avres);
+ if (seg[k]->mdtHitsOnTrack()>3) {
+
+ if (smoothing.addResidualsFromSegment(
+ *seg[k], true, 7.0*avres)) {
+ counter++;
+ }
+
+ }
+ else {
+ if (smoothing.addResidualsFromSegment(*seg[k], false, 7.0*avres)) {
+ counter++;
+ }
+ }
+
+ }
+
+ // break, do no smoothing if there are not enough segments //
+ if (counter<1000) {
+ cerr << "Class RtCalibrationAnalytic, no smoothing applied due to "
+ << "too small number of reconstructed segments!\n";
+ return getResults();
+ }
+
+// // break, do no smoothing if there are not enough segments //
+// if (counter<250) {
+// cerr << "Class RtCalibrationAnalytic, no smoothing applied due to "
+// << "too small number of reconstructed segments!\n";
+// return getResults();
+// }
+//
+// // set bin content for residual bins //
+// unsigned int bin_content(100);
+// if (counter<500) {
+// bin_content = 50;
+// }
+// if (counter>3000) {
+// bin_content = 200;
+// }
+// if (counter>6000) {
+// bin_content = 400;
+// }
+// if (counter>12000) {
+// bin_content = 800;
+// }
+// if (counter>24000) {
+// bin_content = seg.size()/60;
+// }
+//
+// // smoothing //
+// RtRelationLookUp smooth_rt(smoothing.performSmoothing(*tmp_rt, bin_content,
+// m_fix_min, m_fix_max));
+
+ // 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 = sqrt(0.01*RMS);
+
+ // increase the iterations counter //
+ it++;
+
+ // delete tmp_rt and update it //
+// if (it>1) {
+// delete tmp_rt;
+// }
+ delete tmp_rt;
+ tmp_rt = new RtRelationLookUp(smooth_rt);
+// cout << "1: " << tmp_rt << endl;
+
+//---------------------------------------------------------------------------//
+ }
+//---------------------------------------------------------------------------//
+//---------------------------------------------------------------------------//
+
+// // r-t format conversion //
+// vector<double> rt_params;
+// rt_params.push_back(tmp_rt->tLower());
+// rt_params.push_back(0.01*(tmp_rt->tUpper()-tmp_rt->tLower()));
+// for (double t=tmp_rt->tLower(); t<=tmp_rt->tUpper(); t=t+rt_params[1]) {
+// rt_params.push_back(tmp_rt->radius(t));
+// }
+// RtRelationLookUp *improved_rt = new RtRelationLookUp(rt_params);
+
+// ofstream out1("out1.txt");
+// for (double t=conv_rt->tLower(); t<=conv_rt->tUpper(); t=t+5.0) {
+// out1 << t << "\t" << tmp_rt->radius(t) << "\t"
+// << conv_rt->radius(t)
+// << endl;
+// }
+
+//
+// // m_output = new RtCalibrationOutput(smooth_rt,
+// // new RtFullInfo("RtCalibrationAnalytic", m_iteration,
+// // m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+ delete m_output;
+ m_output = new RtCalibrationOutput(tmp_rt,
+ new RtFullInfo("RtCalibrationAnalytic", m_iteration,
+ m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+
+/////////////////////////////////////////
+// RETURN THE RESULT AFTER CONVERGENCE //
+/////////////////////////////////////////
+
+ return getResults();
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD handleSegment ::
+//::::::::::::::::::::::::::
+
+bool RtCalibrationAnalytic::handleSegment(MuonCalibSegment & seg) {
+
+/////////////////////////////////////////////////
+// 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(fabs(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)
+ MTStraightLine track; // refitted straight line
+ vector<double> d_track; // signed distances of the track from the anode
+ // wires of the tubes
+ vector<double> residual_value; // residuals
+ 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
+ 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 = sqrt(av_res/static_cast<double>(seg.mdtHitsOnTrack()));
+ chi2_scale_factor = 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*sqrt(av_res*av_res+
+ m_rt_accuracy*m_rt_accuracy));
+
+// loop over the multilayers //
+
+//-----------------------------------------------------------------------------
+ for (unsigned k=0; k<1+static_cast<unsigned int>(m_split_into_ml);
+ k++) {
+//-----------------------------------------------------------------------------
+
+ if (nb_hits_in_ml[k]<3) {
+ continue;
+ }
+
+// refit the segments within the multilayers //
+ if (!m_tracker.fit(seg, hit_selection[k])) {
+ continue;
+ }
+
+ MTStraightLine track(m_tracker.track());
+ if(std::isnan(track.m_x1()) || std::isnan(track.m_x2()) || std::isnan(track.b_x1()) || std::isnan(track.b_x2()))
+ continue;
+
+
+// check the quality of the fit //
+ if (m_tracker.chi2PerDegreesOfFreedom()>5*chi2_scale_factor) {
+ continue;
+ }
+
+// check whether there are at least three track hits //
+ if (m_control_histograms) {
+ m_nb_segment_hits->Fill(
+ static_cast<double>(m_tracker.numberOfTrackHits()),
+ 1.0);
+ }
+ if (m_tracker.numberOfTrackHits()<3) {
+ continue;
+ }
+
+//reject tracks with silly parameters
+ if (fabs(m_tracker.track().m_x2()) > 8e8) continue;
+
+//for filling into data class
+ if(m_iteration==0)
+ {
+ m_track_slope.Insert(m_tracker.track().m_x2());
+ m_track_position.Insert(m_tracker.track().b_x2());
+ }
+
+// bookkeeping for convergence decision and reliability estimation //
+ m_chi2 = m_chi2+m_tracker.chi2PerDegreesOfFreedom();
+ m_nb_segments_used = m_nb_segments_used+1;
+
+// debug display //
+// display_segment(&seg, display);
+// display << "MESSAGE CONTINUE\n";
+
+// fill the autocalibration objects //
+ // auxiliary variables //
+ track = m_tracker.track(); // refitted track
+ d_track = vector<double>(m_tracker.numberOfTrackHits());
+ residual_value = vector<double>(m_tracker.numberOfTrackHits());
+ w = vector<MTStraightLine>(m_tracker.numberOfTrackHits());
+ G = CLHEP::HepVector(m_tracker.numberOfTrackHits());
+ zeta = vector<double>(m_tracker.numberOfTrackHits());
+
+ // base function values //
+ for (unsigned int l=0; l<m_order; l++) {
+ m_U[l] = CLHEP::HepVector(m_tracker.numberOfTrackHits());
+ for (unsigned int m=0; m<m_tracker.numberOfTrackHits(); m++) {
+ x = (m_tracker.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<m_tracker.numberOfTrackHits(); l++) {
+ w[l] = MTStraightLine(Amg::Vector3D(0.0, (m_tracker.trackHits(
+ )[l]->localPosition()).y(),
+ (m_tracker.trackHits(
+ )[l]->localPosition()).z()),
+ xhat, null, null);
+ d_track[l] = track.signDistFrom(w[l]);
+ residual_value[l] = m_tracker.trackHits()[l]->driftRadius()-
+ fabs(d_track[l]);
+ if (m_control_histograms) {
+ m_residuals->Fill(fabs(d_track[l]), residual_value[l],
+ 1.0);
+ }
+ }
+
+ // loop over all track hits //
+ for (unsigned int l=0; l<m_tracker.numberOfTrackHits(); l++) {
+
+ // analytic calculation of G //
+ for (unsigned int m=0; m<m_tracker.numberOfTrackHits(); m++) {
+ zeta[m] = sqrt(1.0+std::pow(track.m_x2(), 2))*
+ (w[m].positionVector()).z()
+ -track.m_x2()*d_track[m];
+ }
+ for (unsigned int m=0; m<m_tracker.numberOfTrackHits(); m++) {
+ double sum1(0.0), sum2(0.0), sum3(0.0), sum4(0.0);
+ for (unsigned int ll=0; ll <m_tracker.numberOfTrackHits();
+ ll++) {
+ sum1 = sum1+
+ (zeta[l]-zeta[ll])*(zeta[ll]-zeta[m])/
+ (m_tracker.trackHits()[m]->sigma2DriftRadius()
+ *m_tracker.trackHits()[ll]->sigma2DriftRadius()
+ );
+ sum2 = sum2+zeta[ll]/
+ m_tracker.trackHits()[ll]->sigma2DriftRadius();
+ sum3 = sum3+1.0/
+ m_tracker.trackHits()[ll]->sigma2DriftRadius();
+ sum4 = sum4+std::pow(zeta[ll]/
+ m_tracker.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 m_A_tmp(m_A);
+ // autocalibration objects //
+ for (unsigned int p=0; p<m_order; p++) {
+ for (unsigned int pp=p; pp<m_order; pp++) {
+ m_A_tmp[p][pp] = m_A[p][pp]+
+ (dot(G, m_U[p])*dot(G, m_U[pp]))/
+ m_tracker.trackHits(
+ )[l]->sigma2DriftRadius();
+ if(std::isnan(m_A_tmp[p][pp])) return true;
+ }
+ }
+ m_A = m_A_tmp;
+ CLHEP::HepVector m_b_tmp(m_b);
+ for (unsigned int p=0; p<m_order; p++) {
+ m_b_tmp[p] = m_b[p]-residual_value[l]*dot(G, m_U[p])/
+ m_tracker.trackHits(
+ )[l]->sigma2DriftRadius();
+ if(std::isnan(m_b_tmp[p])) return true;
+ }
+ m_b=m_b_tmp;
+ }
+
+//-----------------------------------------------------------------------------
+ }
+//-----------------------------------------------------------------------------
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD setInput ::
+//:::::::::::::::::::::
+
+void RtCalibrationAnalytic::setInput(const IMdtCalibrationOutput * rt_input) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ const RtCalibrationOutput *input(
+ dynamic_cast<const RtCalibrationOutput *>
+ (rt_input));
+
+////////////////////////////////////////////
+// CHECK IF THE OUTPUT CLASS IS SUPPORTED //
+////////////////////////////////////////////
+
+ if (input==0) {
+ cerr << endl
+ << "Class RtCalibrationAnalytic, "
+ << "method setInput: ERROR!\n"
+ << "Calibration input class not supported.\n";
+ exit(1);
+ }
+
+/////////////////////////////////////////////////////////////////
+// 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 //
+ const RtChebyshev *rt_Chebyshev(
+ dynamic_cast<const RtChebyshev *>(m_rt));
+ const RtRelationLookUp *rt_LookUp(
+ dynamic_cast<const RtRelationLookUp *>(m_rt));
+
+ if (rt_Chebyshev==0 && rt_LookUp==0) {
+ cerr << endl
+ << "Class RtCalibrationAnalytic, "
+ << "method setInput: ERROR!\n"
+ << "r-t class not supported.\n";
+ exit(1);
+ }
+
+ // RtChebyshev //
+ if (rt_Chebyshev!=0) {
+ 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!=0) {
+ 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;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::
+//:: METHOD analyse ::
+//::::::::::::::::::::
+
+bool RtCalibrationAnalytic::analyse(void) {
+
+ if(m_tfile!=NULL)
+ 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
+ const RtChebyshev *rt_Chebyshev(
+ dynamic_cast<const RtChebyshev *>(m_rt));
+ const RtRelationLookUp *rt_LookUp(
+ dynamic_cast<const RtRelationLookUp *>(m_rt));
+ double r_corr; // radial correction
+ vector<double> rt_param(m_rt->nPar()); // parameters for the new r-t
+ double x; // reduced time
+ RtParabolicExtrapolation rt_extrapolator; // r-t extrapolator
+ RtFromPoints rt_from_points; // r-t from points
+
+////////////////////////////////////////
+// SOLVE THE AUTOCALIBRATION EQUATION //
+////////////////////////////////////////
+
+ m_alpha = m_A.inverse(ifail)*m_b;
+ if (ifail!=0) {
+ cerr << endl
+ << "Class RtCalibrationAnalytic, method analyse: ERROR!"
+ << endl
+ << "Could not solve the autocalibration equation!"
+ << endl;
+ return false;
+ }
+
+
+
+////////////////////////////////////////
+// CALCULATE THE NEW r-t RELATIONSHIP //
+////////////////////////////////////////
+
+// input r-t is of type RtChebyshev //
+ if (rt_Chebyshev!=0) {
+
+ // set the number of points //
+ if (rt_Chebyshev->numberOfRtParameters()>30) {
+ nb_points = rt_Chebyshev->numberOfRtParameters();
+ }
+
+ // r step size //
+ step = m_r_max/static_cast<double>(nb_points);
+
+ // sample points and Chebyshev fitter //
+ vector<SamplePoint> x_r(nb_points+1);
+ BaseFunctionFitter fitter(rt_Chebyshev->numberOfRtParameters());
+ 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->get_reduced_time(
+ 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) {
+ TGraphErrors *gre= new 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->numberOfRtParameters();
+ k++) {
+ rt_param[k+2] = fitter.coefficients()[k];
+ }
+
+ if (m_rt_new==0) {
+ delete m_rt_new;
+ }
+ m_rt_new = new RtChebyshev(rt_param);
+
+ // parabolic extrapolation // commented, as extr. to end-point is done after convergence
+/* if (m_do_parabolic_extrapolation) {
+ RtRelationLookUp tmprt(performParabolicExtrapolation(false, true,
+ *m_rt_new));
+ if (m_rt_new!=0) {
+ delete m_rt_new;
+ }
+
+ vector<SamplePoint> t_r(nb_points+1);
+ for (unsigned int k=0; k<nb_points+1; k++) {
+ t_r[k].set_x1(t_from_r(k*step));
+ t_r[k].set_x2(tmprt.radius(k*step));
+ t_r[k].set_error(1.0);
+ }
+ m_rt_new = new RtChebyshev(rt_from_points.getRtChebyshev(t_r,
+ rt_Chebyshev->numberOfRtParameters()));
+ }*/
+}
+
+// input-rt is of type RtRelationLookUp //
+ if (rt_LookUp!=0) {
+
+ 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;
+ }
+
+ TGraph *gr(NULL);
+ if(m_control_histograms)
+ gr=new 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];
+ }
+ }
+ }
+
+ if (m_rt_new==0) {
+ delete m_rt_new;
+ }
+ m_rt_new = new 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());
+
+ // parabolic extrapolation //
+// if (m_do_parabolic_extrapolation) {
+// RtRelationLookUp tmprt(performParabolicExtrapolation(false, true,
+// *m_rt_new));
+// if (m_rt_new!=0) {
+// delete m_rt_new;
+// }
+// m_rt_new = new RtRelationLookUp(tmprt);
+// }
+//
+ }
+
+/////////////////////////////////////////////////////////
+// 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 (fabs(r_corr) > r_corr_max){
+ r_corr_max = fabs(r_corr);
+ }
+ m_rt_accuracy = m_rt_accuracy+r_corr*r_corr;
+ }
+ m_rt_accuracy = 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);
+ // cout << "------------------ ITERATION: "<< m_iteration <<" ---------------------- " << endl;
+ // cout << " CHI2:\t" << m_chi2 <<endl;
+ // cout << " CORDUNC.DIFF:\t" << fabs(m_rt_accuracy_diff) << endl;
+ // cout << " Accuracy Guess:\t" << fabs(m_rt_accuracy) << endl;
+ // if ((m_chi2<=m_chi2_previous || fabs(m_rt_accuracy_diff)>0.001) &&
+ // m_iteration<m_max_it) {
+ if ( ((m_chi2<m_chi2_previous && fabs(m_chi2-m_chi2_previous)>0.01)
+ || fabs(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;
+
+ if (m_output!=0) {
+ delete m_output;
+ }
+ m_output = new RtCalibrationOutput(m_rt_new,
+ new RtFullInfo("RtCalibrationAnalytic", m_iteration,
+ m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD converged ::
+//::::::::::::::::::::::
+
+bool RtCalibrationAnalytic::converged(void) const {
+
+ return (m_status>0);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD getResults ::
+//:::::::::::::::::::::::
+
+const IMdtCalibrationOutput * RtCalibrationAnalytic::getResults(void) const {
+
+ return static_cast<const IMdtCalibrationOutput *>(m_output);
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD performParabolicExtrapolation ::
+//::::::::::::::::::::::::::::::::::::::::::
+
+RtRelationLookUp * RtCalibrationAnalytic::performParabolicExtrapolation(
+ const bool & min,
+ const bool & max,
+ const IRtRelation & in_rt) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ RtParabolicExtrapolation rt_extrapolator; // r-t extrapolator
+ RtRelationLookUp *rt_low(0), *rt_high(0); // pointers to the r-t
+ // relationships after
+ // extrapolation
+ 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 = new 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 = new 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_low = new RtRelationLookUp(
+ rt_extrapolator.getRtWithParabolicExtrapolation(*rt_high,
+ 1.0,
+ 3.0,
+ 0.0,
+ add_fit_point));
+ } else {
+ rt_low = new RtRelationLookUp(
+ rt_extrapolator.getRtWithParabolicExtrapolation(in_rt,
+ 1.0,
+ 3.0,
+ 0.0,
+ add_fit_point));
+ }
+ }
+
+////////////////////
+// RETURN RESULTS //
+////////////////////
+
+ if (min && max) {
+ if(rt_high) delete rt_high;
+ return rt_low;
+ }
+ if (min) {
+ if(rt_high) delete rt_high;
+ return rt_low;
+ }
+ if (rt_low) delete rt_low;
+ return rt_high;
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx
new file mode 100644
index 0000000000000..9a2b92d5868f2
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationCurved.cxx
@@ -0,0 +1,1699 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+ //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 10.05.2008, AUTHOR: OLIVER KORTNER
+// Modified: 01.03.2009 by O. Kortner, parabolic extrapolation added.
+// 15.03.2009 by O. Kortner, smoothing added.
+// 05.02.2010, JOERG v. LOEBEN
+// major changes: - parabolic extrapolation to r(t_max) now done after the
+// after the method converged and - if turned on - after the
+// the residual smoothing.
+// - prevent extrapolation to r_max to become too large.
+// - changed convergence criterion of the residual smoothing
+// from 1 mum to 2 mum (tuned on cosmics)
+// - reduced the roadwidth of the curved segment refit for
+// smoothing from 7->5 times the average
+// (tuned on cosmics) tube resolution
+// - change in convergence criterion to prevent
+// iterating with only minor changes (tuned on cosmics)
+// - increase of spacial precision in method t_from_r to 1mum
+// (important to reduce oscillations in r(t) relationship)
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS RtCalibrationCurved ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibRt/RtCalibrationCurved.h"
+#include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
+#include "MdtCalibRt/RtCalibrationOutput.h"
+#include "MdtCalibData/RtChebyshev.h"
+#include "MdtCalibData/RtRelationLookUp.h"
+#include "MuonCalibMath/BaseFunctionFitter.h"
+#include "MuonCalibMath/ChebyshevPolynomial.h"
+#include "MuonCalibMath/LegendrePolynomial.h"
+#include "MuonCalibMath/PolygonBase.h"
+#include "CLHEP/Matrix/Matrix.h"
+#include "MdtCalibRt/RtParabolicExtrapolation.h"
+#include "MdtCalibData/RtFromPoints.h"
+#include "MdtCalibRt/AdaptiveResidualSmoothing.h"
+
+#include "MdtCalibInterfaces/IMdtCalibrationOutput.h"
+#include "MdtCalibData/IRtRelation.h"
+#include "MdtCalibData/RtRelationLookUp.h"
+#include "MdtCalibData/RtScaleFunction.h"
+#include "MdtCalibRt/RtCalibrationOutput.h"
+#include "MdtCalibFitters/CurvedPatRec.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "MuonCalibMath/BaseFunction.h"
+
+#include "MultilayerRtDifference.h"
+
+
+#include "sstream"
+#include "TF1.h"
+#include "TProfile.h"
+
+//::::::::::::::::::::::::
+//:: NAMESPACE SETTINGS ::
+//::::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+inline Double_t RtCalibrationCurved_polfun(Double_t *x, Double_t *par)
+ {
+ return par[0] * RtScalePolynomial(x[0]);
+ }
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: DEFAULT CONSTRUCTOR ::
+//:::::::::::::::::::::::::
+
+RtCalibrationCurved::RtCalibrationCurved(std::string name) :
+ IMdtCalibration(name), m_rt_accuracy_previous(0.0){
+
+ init(0.5*CLHEP::mm, 1, 15, true, false, 15, false, false, false);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: CONSTRUCTOR ::
+//:::::::::::::::::
+
+RtCalibrationCurved::RtCalibrationCurved(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,
+ bool do_smoothing,
+ bool do_multilayer_rt_scale)
+ : 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);
+
+}
+
+
+//*****************************************************************************
+
+//::::::::::::::::
+//:: DESTRUCTOR ::
+//::::::::::::::::
+
+RtCalibrationCurved::~RtCalibrationCurved(void) {
+
+ if (m_tracker!=0) {
+ delete m_tracker;
+ }
+
+ if (m_base_function!=0) {
+ delete m_base_function;
+ }
+ delete m_multilayer_rt_difference;
+ if (m_tfile!=0) {
+ m_tfile->Write();
+ delete m_cut_evolution;
+ delete m_nb_segment_hits;
+ delete m_pull_initial;
+ delete m_pull_final;
+ delete m_residuals_initial;
+ delete m_residuals_final;
+ delete m_tfile;
+ }
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD reliability ::
+//::::::::::::::::::::::::
+
+double RtCalibrationCurved::reliability(void) const {
+
+ return m_status;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::
+//:: METHOD estimatedRtAccuracy ::
+//::::::::::::::::::::::::::::::::
+
+double RtCalibrationCurved::estimatedRtAccuracy(void) const {
+
+ return m_rt_accuracy;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::
+//:: METHOD numberOfSegments ::
+//:::::::::::::::::::::::::::::
+
+int RtCalibrationCurved::numberOfSegments(void) const {
+
+ return m_nb_segments;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::
+//:: METHOD numberOfSegmentsUsed ::
+//:::::::::::::::::::::::::::::::::
+
+int RtCalibrationCurved::numberOfSegmentsUsed(void) const {
+
+ return m_nb_segments_used;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD iteration ::
+//::::::::::::::::::::::
+
+int RtCalibrationCurved::iteration(void) const {
+
+ return m_iteration;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD smoothing ::
+//::::::::::::::::::::::
+
+bool RtCalibrationCurved::smoothing(void) const {
+
+ return m_do_smoothing;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::
+//:: METHOD setEstimateRtAccuracy ::
+//::::::::::::::::::::::::::::::::::
+
+void RtCalibrationCurved::setEstimateRtAccuracy(const double & acc) {
+
+ m_rt_accuracy = fabs(acc);
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::::::
+//:: METHOD switch_on_control_histograms ::
+//:::::::::::::::::::::::::::::::::::::::::
+
+void RtCalibrationCurved::switch_on_control_histograms(
+ const std::string & file_name) {
+
+/////////////////////////////////////////////
+// CREATE THE ROOT FILE AND THE HISTOGRAMS //
+/////////////////////////////////////////////
+
+ m_control_histograms = true;
+
+ m_tfile = new TFile(file_name.c_str(), "RECREATE");
+
+ m_cut_evolution = new TH1F("m_cut_evolution",
+ "CUT EVOLUTION",
+ 11, -0.5, 10.5);
+
+ m_nb_segment_hits = new TH1F("m_nb_segment_hits",
+ "NUMBER OF HITS ON THE REFITTED SEGMENTS",
+ 11, -0.5, 10.5);
+
+ m_pull_initial = new TH1F("m_pull_initial",
+ "INITIAL PULL DISTRIBUTION",
+ 200, -5.05, 5.05);
+// m_pull_final = new TH1F("m_pull_final",
+// "FINAL PULL DISTRIBUTION",
+// 200, -5.05, 5.05);
+
+ m_residuals_initial = new TH2F("m_residuals_initial",
+ "INITIAL OF THE REFITTED SEGMENTS",
+ 100, -0.5, 15.0, 300, -1.5, 1.5);
+ m_residuals_final = new TH2F("m_residuals_final",
+ "FINAL OF THE REFITTED SEGMENTS",
+ 100, -0.5, 15.0, 300, -1.5, 1.5);
+
+ delete m_multilayer_rt_difference;
+ m_multilayer_rt_difference=new MultilayerRtDifference(10000, m_tfile);
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD switch_off_control_histograms ::
+//::::::::::::::::::::::::::::::::::::::::::
+
+void RtCalibrationCurved::switch_off_control_histograms(void) {
+
+ m_control_histograms = false;
+ if (m_tfile!=0) {
+ m_tfile->Write();
+ delete m_multilayer_rt_difference;
+ m_multilayer_rt_difference=new MultilayerRtDifference(10000);
+ }
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD forceMonotony ::
+//::::::::::::::::::::::::::
+
+void RtCalibrationCurved::forceMonotony(void) {
+
+ m_force_monotony = true;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::
+//:: METHOD doNotForceMonotony ::
+//:::::::::::::::::::::::::::::::
+
+void RtCalibrationCurved::doNotForceMonotony(void) {
+
+ m_force_monotony = false;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::
+//:: METHOD doParabolicExtrapolation ::
+//:::::::::::::::::::::::::::::::::::::
+
+void RtCalibrationCurved::doParabolicExtrapolation(void) {
+
+ m_do_parabolic_extrapolation = true;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD doSmoothing ::
+//::::::::::::::::::::::::::
+
+void RtCalibrationCurved::doSmoothing(void) {
+
+ m_do_smoothing = true;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD noSmoothing ::
+//::::::::::::::::::::::::::
+
+void RtCalibrationCurved::noSmoothing(void) {
+
+ m_do_smoothing = false;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::
+//:: METHOD analyseSegments ::
+//::::::::::::::::::::::::::::
+
+const IMdtCalibrationOutput * RtCalibrationCurved::analyseSegments(
+ const std::vector<MuonCalibSegment*> & seg) {
+
+ const IRtRelation *tmp_rt(0);
+
+//////////////////////////
+// AUTOCALIBRATION LOOP //
+//////////////////////////
+
+ while (!converged()) {
+
+ for (unsigned int k=0; k<seg.size(); k++) {
+ handleSegment(*seg[k]);
+ }
+ if(!analyse(seg)) {
+ for(unsigned int i=0; i<seg.size(); i++) {
+ std::cout << i << " " << seg[i]->direction() <<" "
+ << seg[i]->position() << std::endl;
+ }
+ return NULL;
+ }
+
+ const RtCalibrationOutput* rtOut =m_output;
+
+ if (!converged()) {
+ setInput(rtOut);
+ }
+ tmp_rt = rtOut->rt();
+ }
+
+// parabolic extrapolations for small radii //
+ if (m_do_parabolic_extrapolation) {
+ RtRelationLookUp *tmprt =
+ performParabolicExtrapolation(true, true, *tmp_rt);
+ delete m_output;
+ m_output = new RtCalibrationOutput(tmprt,
+ new RtFullInfo("RtCalibrationCurved", m_iteration,
+ m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+ delete tmp_rt;
+ tmp_rt = tmprt;
+ }
+
+//////////////////////////////////////////////
+// SMOOTHING AFTER CONVERGENCE IF REQUESTED //
+//////////////////////////////////////////////
+
+ if (!m_do_smoothing) {
+// final residuals //
+ double r, d;
+ for (unsigned int k=0; k<seg.size(); k++) {
+ for (unsigned int l=0; l<seg[k]->hitsOnTrack(); l++) {
+ r = fabs((seg[k]->mdtHOT())[l]->driftRadius());
+ d = fabs((seg[k]->mdtHOT())[l]->signedDistanceToTrack());
+ if(m_residuals_final != NULL)
+ m_residuals_final->Fill(d, r-d, 1.0);
+ }
+ }
+ 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);
+ unsigned int counter2(0);
+
+ // overwrite drift radii and calculate the average resolution //
+ for (unsigned int k=0; k<seg.size(); k++) {
+ if (seg[k]->mdtHitsOnTrack()<4) {
+ continue;
+ }
+ counter2++;
+ double avres(0.0);
+ for (unsigned int h=0; h<seg[k]->mdtHitsOnTrack(); h++) {
+ seg[k]->mdtHOT()[h]->setDriftRadius(tmp_rt->radius(
+ seg[k]->mdtHOT()[h]->driftTime()),
+ seg[k]->mdtHOT()[h]->sigmaDriftRadius());
+ if (seg[k]->mdtHOT()[h]->sigmaDriftRadius()<0.5*m_r_max) {
+ avres = avres+seg[k]->mdtHOT()[h]->sigma2DriftRadius();
+ } else {
+ avres = avres+0.01;
+ }
+ }
+ avres = avres/static_cast<double>(seg[k]->mdtHitsOnTrack());
+ avres = sqrt(avres);
+ if (smoothing.addResidualsFromSegment(*seg[k], true, 5.0*avres)) {
+ counter++;
+ }
+ }
+
+ // break, do no smoothing if there are not enough segments //
+ if (counter<1000) {
+ cerr << "Class RtCalibrationCurved, no smoothing applied due to "
+ << "too small number of reconstructed segments!\n";
+// final residuals //
+ double r, d;
+ for (unsigned int k=0; k<seg.size(); k++) {
+ for (unsigned int l=0; l<seg[k]->hitsOnTrack(); l++) {
+ r = fabs((seg[k]->mdtHOT())[l]->driftRadius());
+ d = fabs((seg[k]->mdtHOT())[l]->signedDistanceToTrack());
+ if(m_residuals_final != NULL)
+ m_residuals_final->Fill(d, r-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 = sqrt(0.01*RMS);
+
+ // increase the iterations counter //
+ it++;
+
+ // delete tmp_rt and update it //
+ delete tmp_rt;
+ tmp_rt = new RtRelationLookUp(smooth_rt);
+
+//---------------------------------------------------------------------------//
+ }
+//---------------------------------------------------------------------------//
+//---------------------------------------------------------------------------//
+
+ delete m_output;
+ m_output = new RtCalibrationOutput(tmp_rt,
+ new RtFullInfo("RtCalibrationCurved", m_iteration,
+ m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+
+/////////////////////
+// FINAL RESIDUALS //
+/////////////////////
+
+ double r, d;
+ for (unsigned int k=0; k<seg.size(); k++) {
+ for (unsigned int l=0; l<seg[k]->hitsOnTrack(); l++) {
+ r = fabs((seg[k]->mdtHOT())[l]->driftRadius());
+ d = fabs((seg[k]->mdtHOT())[l]->signedDistanceToTrack());
+ if(m_residuals_final != NULL)
+ m_residuals_final->Fill(d, r-d, 1.0);
+ }
+ }
+
+/////////////////////////////////////////
+// RETURN THE RESULT AFTER CONVERGENCE //
+/////////////////////////////////////////
+
+ return getResults();
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD handleSegment ::
+//::::::::::::::::::::::::::
+
+bool RtCalibrationCurved::handleSegment(MuonCalibSegment & seg) {
+
+/////////////////////////////////////////////////
+// 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 (fabs(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)
+ vector<double> d_track; // signed distances of the track from the anode
+ // wires of the tubes
+ 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 //
+ 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 //
+///////////////////////////////////////
+
+// 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 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 //
+// std::cout<<"XXxxXX "<<m_rt->GetTmaxDiff()<<std::endl;
+ 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 = sqrt(av_res/static_cast<double>(seg.mdtHitsOnTrack()));
+ chi2_scale_factor = 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*sqrt(av_res*av_res+m_rt_accuracy*m_rt_accuracy));
+
+// check whether there are enough hits in the chambers //
+ if (nb_hits_in_ml<4) {
+ return true;
+ }
+
+// refit the segments //
+ if (!m_tracker->fit(seg, hit_selection)) {
+ return true;
+ }
+
+ CurvedLine track(m_tracker->curvedTrack());
+ if(std::isnan(m_tracker->chi2())) {
+ return true;
+ }
+
+// check the quality of the fit //
+ if (m_tracker->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(
+ static_cast<double>(m_tracker->numberOfTrackHits()), 1.0);
+ }
+ if (m_tracker->numberOfTrackHits()<4) {
+ return true;
+ }
+
+// display_segment(&seg, display&(m_tracker->curvedTrack()));
+
+//reject tracks with silly parameters
+ if (fabs(m_tracker->curvedTrack().getTangent(
+ seg.mdtHOT()[0]->localPosition().z()).m_x2())>8.0e8) {
+ return true;
+ }
+
+// bookkeeping for convergence decision and reliability estimation //
+ m_chi2 = m_chi2+m_tracker->chi2PerDegreesOfFreedom();
+ m_nb_segments_used = m_nb_segments_used+1;
+
+// debug display //
+// display_segment(&seg, display);
+// display << "MESSAGE CONTINUE\n";
+
+// fill the autocalibration objects //
+
+
+ // track coeffient matrix and its inverse //
+ F = vector<CLHEP::HepVector>(m_tracker->numberOfTrackHits());
+ for (unsigned int h=0; h<m_tracker->numberOfTrackHits(); h++) {
+
+ const MdtCalibHitBase &hb=*(m_tracker->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++) {
+ (F[h])[p] = m_Legendre->value(p, (m_tracker->trackHits(
+ )[h]->localPosition()).z())
+ /sqrt(1.0+std::pow(track.getTangent(
+ (m_tracker->trackHits(
+ )[h]->localPosition()).z()).m_x2(), 2));
+ }
+ }
+ 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<m_tracker->numberOfTrackHits(); h++) {
+ m_M_track[p][q] = m_M_track[p][q]+(F[h])[p]*(F[h])[q]/
+ (m_tracker->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) {
+ cerr << endl
+ << "Class RtCalibrationCurved, method handleSegment: WARNING!"
+ << endl
+ << "Could not invert track matrix!\n";
+ return true;
+ }
+
+// Jacobian matrix for the residuals //
+ // track distances, residuals, corrections //
+ d_track = vector<double>(m_tracker->numberOfTrackHits());
+ w = vector<MTStraightLine>(m_tracker->numberOfTrackHits());
+ residual_value = CLHEP::HepVector(m_tracker->numberOfTrackHits());
+ weighted_residual = CLHEP::HepVector(m_tracker->numberOfTrackHits());
+ for (unsigned int l=0; l<m_order; l++) {
+ m_U[l] = CLHEP::HepVector(m_tracker->numberOfTrackHits());
+ m_U_weighted[l] = CLHEP::HepVector(m_tracker->numberOfTrackHits());
+ }
+ for (unsigned int h=0; h<m_tracker->numberOfTrackHits(); h++) {
+ w[h] = MTStraightLine(Amg::Vector3D(0.0, (m_tracker->trackHits(
+ )[h]->localPosition()).y(),
+ (m_tracker->trackHits(
+ )[h]->localPosition()).z()), xhat, null, null);
+ d_track[h] = track.getTangent((m_tracker->trackHits(
+ )[h]->localPosition()).z()).signDistFrom(w[h]);
+ residual_value[h] = (m_tracker->trackHits()[h])->driftRadius()-
+ fabs(d_track[h]);
+ if (m_control_histograms) {
+ if (m_iteration==0) {
+ m_residuals_initial->Fill(fabs(d_track[h]), residual_value[h],
+ 1.0);
+ }
+ }
+ for (unsigned int l=0; l<m_order; l++) {
+ x = (m_tracker->trackHits()[h]->driftRadius()
+ -0.5*m_r_max)/(0.5*m_r_max);
+ (m_U[l])[h] = m_base_function->value(l, x);
+ }
+ }
+
+
+// test track matrix //
+// static ofstream outfile("tmp.txt");
+// CLHEP::HepVector my_b(3);
+// for (int p=0; p<my_b.num_row(); p++) {
+// my_b[p] = 0.0;
+// for (unsigned int h=0; h<m_tracker->numberOfTrackHits(); h++) {
+// // my_b[p] = my_b[p]+m_Legendre->value(p, (m_tracker->trackHits(
+// // )[h]->localPosition()).z())*
+// // (
+// // +(2*d_track[h]>=0-1)*m_tracker->trackHits()[h]->driftRadius()
+// // +(m_tracker->trackHits()[h]->localPosition()).y()
+// // /sqrt(1.0+std::pow((track.getTangent(
+// // (m_tracker->trackHits(
+// // )[h]->localPosition()).z())).m_x2(), 2))
+// // )
+// // /(
+// // (m_tracker->trackHits()[h])->sigma2DriftRadius()
+// // *sqrt(1.0+std::pow((track.getTangent(
+// // (m_tracker->trackHits(
+// // )[h]->localPosition()).z())).m_x2(), 2))
+// // );
+// my_b[p] = my_b[p]+((m_tracker->trackHits()[h]->localPosition()).y()
+// +(2*(d_track[h]>=0)-1)
+// *sqrt(1.0+std::pow((track.getTangent(
+// (m_tracker->trackHits(
+// )[h]->localPosition()).z())).m_x2() ,2))
+// *m_tracker->trackHits()[h]->driftRadius())
+// *m_Legendre->value(p, (m_tracker->trackHits(
+// )[h]->localPosition()).z())
+// /((m_tracker->trackHits()[h])->sigma2DriftRadius()
+// *(1.0+std::pow((track.getTangent(
+// (m_tracker->trackHits(
+// )[h]->localPosition()).z())).m_x2() ,2)));
+// }
+// }
+// CLHEP::HepVector my_alpha(3);
+// my_alpha = m_M_track_inverse*my_b;
+// for (unsigned int h=0; h<m_tracker->numberOfTrackHits(); h++) {
+// double ypos(0.0);
+// for (unsigned int p=0; p<3; p++) {
+// ypos = ypos+my_alpha[p]*m_Legendre->value(p, (m_tracker->trackHits(
+// )[h]->localPosition()).z());
+// }
+// outfile << (m_tracker->trackHits()[h]->localPosition()).z() << "\t"
+// << m_tracker->curvedTrack().getPointOnLine(
+// (m_tracker->trackHits()[h]->localPosition()).z()).y()
+// << "\t" << ypos << endl;
+// }
+
+ // Jacobian //
+ D = CLHEP::HepMatrix(m_tracker->numberOfTrackHits(), m_tracker->numberOfTrackHits());
+ for (unsigned int h=0; h<m_tracker->numberOfTrackHits(); h++) {
+ for (unsigned int hp=0; hp<m_tracker->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])/
+ (m_tracker->trackHits()[hp])->sigma2DriftRadius();
+ }
+ }
+
+// autocalibration objects //
+ // errors of the residuals //
+ vector<double> sigma_residual(m_tracker->numberOfTrackHits(), 0.0);
+ for (unsigned int h=0; h<m_tracker->numberOfTrackHits(); h++) {
+ for (unsigned int hp=0; hp<m_tracker->numberOfTrackHits(); hp++) {
+ sigma_residual[h] = sigma_residual[h]+
+ std::pow(D[h][hp]*(m_tracker->trackHits()[hp]
+ )->sigmaDriftRadius(), 2);
+ }
+ sigma_residual[h] = sqrt(sigma_residual[h]);
+ if (sigma_residual[h]<av_res/
+ static_cast<double>(m_tracker->numberOfTrackHits())) {
+ sigma_residual[h] = av_res/sqrt(m_tracker->numberOfTrackHits());
+ }
+ }
+
+
+ // weight residuals and correction functions //
+ for (unsigned int h=0; h<m_tracker->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 m_A_tmp(m_A);
+
+ // autocalibration objects //
+ for (unsigned int p=0; p<m_order; p++) {
+ for (unsigned int pp=p; pp<m_order; pp++) {
+ m_A_tmp[p][pp] = m_A[p][pp]+dot(D*m_U_weighted[p], D*m_U_weighted[pp]);
+ if (std::isnan(m_A_tmp[p][pp])) {
+ return true;
+ }
+ }
+ }
+
+ CLHEP::HepVector m_b_tmp(m_b);
+ for (unsigned int p=0; p<m_order; p++) {
+ m_b_tmp[p] = m_b[p]+dot(D*m_U_weighted[p], weighted_residual);
+ if (std::isnan(m_b_tmp[p])) {
+ return true;
+ }
+ }
+
+ m_A = m_A_tmp;
+ m_b=m_b_tmp;
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD setInput ::
+//:::::::::::::::::::::
+
+void RtCalibrationCurved::setInput(const IMdtCalibrationOutput * rt_input) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ const RtCalibrationOutput *input(
+ dynamic_cast<const RtCalibrationOutput *>
+ (rt_input));
+
+////////////////////////////////////////////
+// CHECK IF THE OUTPUT CLASS IS SUPPORTED //
+////////////////////////////////////////////
+
+ if (input==0) {
+ cerr << endl
+ << "Class RtCalibrationCurved, "
+ << "method setInput: ERROR!\n"
+ << "Calibration input class not supported.\n";
+ exit(1);
+ }
+
+/////////////////////////////////////////////////////////////////
+// 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 //
+ const RtChebyshev *rt_Chebyshev(
+ dynamic_cast<const RtChebyshev *>(m_rt));
+ const RtRelationLookUp *rt_LookUp(
+ dynamic_cast<const RtRelationLookUp *>(m_rt));
+
+ if (rt_Chebyshev==0 && rt_LookUp==0) {
+ cerr << endl
+ << "Class RtCalibrationCurved, "
+ << "method setInput: ERROR!\n"
+ << "r-t class not supported.\n";
+ exit(1);
+ }
+
+ // RtChebyshev //
+ if (rt_Chebyshev!=0) {
+ 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!=0) {
+ 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;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::
+//:: METHOD analyse ::
+//::::::::::::::::::::
+
+bool RtCalibrationCurved::analyse(const std::vector<MuonCalibSegment*> & seg) {
+
+///////////////
+// 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
+ const RtChebyshev *rt_Chebyshev(
+ dynamic_cast<const RtChebyshev *>(m_rt));
+ const RtRelationLookUp *rt_LookUp(
+ dynamic_cast<const RtRelationLookUp *>(m_rt));
+ double r_corr; // radial correction
+ vector<double> rt_param(m_rt->nPar()); // parameters for the new r-t
+ double x; // reduced time
+ RtFromPoints rt_from_points; // r-t from points
+
+////////////////////////////////////////
+// SOLVE THE AUTOCALIBRATION EQUATION //
+////////////////////////////////////////
+
+ m_alpha = m_A.inverse(ifail)*m_b;
+ if (ifail!=0) {
+ cerr << endl
+ << "Class RtCalibrationCurved, method analyse: ERROR!"
+ << endl
+ << "Could not solve the autocalibration equation!"
+ << endl;
+ return false;
+ }
+
+////////////////////////////////////////
+// CALCULATE THE NEW r-t RELATIONSHIP //
+////////////////////////////////////////
+
+// input r-t is of type RtChebyshev //
+ if (rt_Chebyshev!=0) {
+
+ // set the number of points //
+ if (rt_Chebyshev->numberOfRtParameters()>30) {
+ nb_points = rt_Chebyshev->numberOfRtParameters();
+ }
+
+ // r step size //
+ step = m_r_max/static_cast<double>(nb_points);
+
+ // sample points and Chebyshev fitter //
+ vector<SamplePoint> x_r(nb_points+1);
+ BaseFunctionFitter fitter(rt_Chebyshev->numberOfRtParameters());
+ 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->get_reduced_time(
+ 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());
+ }
+ }
+ }
+
+ // 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->numberOfRtParameters();
+ k++) {
+ rt_param[k+2] = fitter.coefficients()[k];
+ }
+
+ if (m_rt_new==0) {
+ delete m_rt_new;
+ }
+ m_rt_new = new RtChebyshev(rt_param);
+ m_rt_new->SetTmaxDiff(m_rt->GetTmaxDiff());
+
+ // parabolic extrapolation //
+ /*
+ if (m_do_parabolic_extrapolation) {
+ RtRelationLookUp tmprt(performParabolicExtrapolation(false, true,
+ *m_rt_new));
+// vector<SamplePoint> add_fit_point;
+// if (m_fix_min) {
+// add_fit_point.push_back(
+// SamplePoint(m_rt_new->tLower(), 0.0, 1.0));
+// }
+// RtRelationLookUp tmp_rt(rt_extrapolator.getRtWithParabolicExtrapolation(
+// *m_rt_new,
+// m_rt_new->radius(m_rt_new->tLower())+1.0,
+// m_rt_new->radius(m_rt_new->tLower())+3.0,
+// m_rt_new->radius(m_rt_new->tLower()),
+// add_fit_point));
+// add_fit_point.clear();
+// if (m_fix_max) {
+// add_fit_point.push_back(
+// SamplePoint(m_rt_new->tUpper(), m_r_max, 1.0));
+// }
+// RtRelationLookUp tmprt(rt_extrapolator.getRtWithParabolicExtrapolation(
+// tmp_rt,
+// m_rt_new->radius(m_rt_new->tUpper())-1.0,
+// m_rt_new->radius(m_rt_new->tUpper())-3.0,
+// m_rt_new->radius(m_rt_new->tUpper()),
+// add_fit_point));
+
+ if (m_rt_new!=0) {
+ delete m_rt_new;
+ }
+
+ vector<SamplePoint> t_r(nb_points+1);
+ for (unsigned int k=0; k<nb_points+1; k++) {
+ t_r[k].set_x1(t_from_r(k*step));
+ t_r[k].set_x2(tmprt.radius(k*step));
+ t_r[k].set_error(1.0);
+ }
+ m_rt_new = new RtChebyshev(rt_from_points.getRtChebyshev(t_r,
+ rt_Chebyshev->numberOfRtParameters()));
+ }
+*/
+ }
+
+// input-rt is of type RtRelationLookUp //
+ if (rt_LookUp!=0) {
+
+ 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] - 7.6)/7.6;
+ 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];
+ }
+ }
+ }
+
+ if (m_rt_new==0) {
+ delete m_rt_new;
+ }
+ m_rt_new = new RtRelationLookUp(rt_param);
+ m_rt_new->SetTmaxDiff(m_rt->GetTmaxDiff());
+ // parabolic extrapolation //
+/*
+ if (m_do_parabolic_extrapolation) {
+ RtRelationLookUp tmprt(performParabolicExtrapolation(false, true,
+ *m_rt_new));
+// vector<SamplePoint> add_fit_point;
+// // if (m_fix_min) {
+// // add_fit_point.push_back(
+// // SamplePoint(m_rt_new->tLower(), 0.0, 1.0));
+// // }
+// // // RtRelationLookUp tmp_rt(rt_extrapolator.getRtWithParabolicExtrapolation(
+// // // *m_rt_new,
+// // // m_rt_new->radius(m_rt_new->tLower())+1.0,
+// // // m_rt_new->radius(m_rt_new->tLower())+3.0,
+// // // m_rt_new->radius(m_rt_new->tLower()),
+// // // add_fit_point));
+// // RtRelationLookUp tmp_rt(rt_extrapolator.getRtWithParabolicExtrapolation(
+// // *m_rt_new,
+// // 1.2,
+// // 3.0,
+// // 0.0,
+// // add_fit_point));
+// // add_fit_point.clear();
+// if (m_fix_max) {
+// add_fit_point.push_back(
+// SamplePoint(m_rt_new->tUpper(), m_r_max, 1.0));
+// }
+// RtRelationLookUp tmprt(rt_extrapolator.getRtWithParabolicExtrapolation(
+// *m_rt_new,
+// m_r_max-1.0,
+// m_r_max-3.0,
+// m_r_max,
+// add_fit_point));
+
+ if (m_rt_new!=0) {
+ delete m_rt_new;
+ }
+ m_rt_new = new RtRelationLookUp(tmprt);
+ }
+*/
+ }
+
+/////////////////////////////////////////////////////////
+// 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 (fabs(r_corr) > r_corr_max){
+ r_corr_max = fabs(r_corr);
+ }
+ m_rt_accuracy = m_rt_accuracy+r_corr*r_corr;
+ }
+ m_rt_accuracy = 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 || fabs(m_chi2-m_chi2_previous)>0.01)
+ || (fabs(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, &seg);
+ else
+ {
+// std::cout<<"No update in first iteration!"<<std::endl;
+ m_multilayer_rt_difference->DoFit();
+ }
+
+//////////////////////////////////////////////////
+// STORE THE RESULTS IN THE RtCalibrationOutput //
+//////////////////////////////////////////////////
+
+ m_iteration = m_iteration+1;
+
+ if (m_output!=0) {
+ delete m_output;
+ }
+ m_output = new RtCalibrationOutput(m_rt_new,
+ new RtFullInfo("RtCalibrationCurved", m_iteration,
+ m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+
+
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD converged ::
+//::::::::::::::::::::::
+
+bool RtCalibrationCurved::converged(void) const {
+
+ return (m_status>0);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD getResults ::
+//:::::::::::::::::::::::
+
+const IMdtCalibrationOutput * RtCalibrationCurved::getResults(void) const {
+
+ return static_cast<const IMdtCalibrationOutput *>(m_output);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD 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) {
+
+/////////////////////////////
+// RESET PRIVATE VARIABLES //
+/////////////////////////////
+ m_rt=NULL;
+ m_r_max = 15.0*CLHEP::mm;
+ m_control_histograms = false;
+ m_tfile = 0;
+ m_cut_evolution = 0;
+ m_nb_segment_hits = 0;
+ m_pull_initial = 0;
+ m_pull_final = 0;
+ m_residuals_initial = 0;
+ m_residuals_final = 0;
+ 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 = fabs(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 = abs(max_it);
+ m_do_multilayer_rt_scale=do_multilayer_rt_scale;
+ m_multilayer_rt_difference = new MultilayerRtDifference(10000);
+
+ m_tracker = new CurvedPatRec;
+
+ if (m_order==0) {
+ cerr << "\n"
+ << "Class RtCalibrationCurved, method init: ERROR!"
+ << "\nOrder of the correction polynomial must be >0!\n";
+ exit(1);
+ }
+
+ 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_rt_new = 0;
+ m_output = 0;
+
+ m_U = vector<CLHEP::HepVector>(m_order);
+ m_U_weighted = 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 = 0;
+ cerr << "\n"
+ << "Class RtCalibrationCurved, method init: "
+ << "ERROR!\n"
+ << "Illegal correction function type!\n";
+ exit(1);
+ }
+ switch(func_type) {
+ case 1:
+ m_base_function = new LegendrePolynomial;
+ break;
+ case 2:
+ m_base_function = new ChebyshevPolynomial;
+ break;
+ case 3:
+ if (m_order<2) {
+ cerr << "\n"
+ << "Class RtCalibrationCurved, "
+ << "method init: ERROR!\n"
+ << "Order must be >2 for polygons! "
+ << "It is set to " << m_order
+ << "by the user.\n";
+ exit(1);
+ }
+ 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 = new 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;
+
+// Legendre polynomials and tracking objects //
+ m_Legendre = Legendre_polynomial::get_Legendre_polynomial();
+ m_M_track = CLHEP::HepSymMatrix(3);
+ m_M_track_inverse = CLHEP::HepSymMatrix(3);
+
+
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD t_from_r ::
+//:::::::::::::::::::::
+
+double RtCalibrationCurved::t_from_r(const double & r) {
+
+///////////////
+// 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 &&
+ fabs(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);
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::
+//:: METHOD display_segment ::
+//::::::::::::::::::::::::::::
+
+void RtCalibrationCurved::display_segment(MuonCalibSegment * segment,
+ std::ofstream & outfile, const CurvedLine * curved_segment) {
+
+///////////////
+// 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 << "NULL " << 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 << "NULL " << 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==0) {
+ MTStraightLine aux_track(segment->position(), segment->direction(),
+ null, null);
+ outfile << "SET PLCI 4\n"
+ << "LINE " << aux_track.m_x2()*(z_min-30.0)+aux_track.b_x2()
+ << " " << z_min-30.0
+ << " " << aux_track.m_x2()*(z_max+30.0)+aux_track.b_x2()
+ << " " << z_max+30.0 << "\n";
+ }
+
+// a curved segment is drawn on demand //
+ if (curved_segment!=0) {
+ 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;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD performParabolicExtrapolation ::
+//::::::::::::::::::::::::::::::::::::::::::
+
+RtRelationLookUp * RtCalibrationCurved::performParabolicExtrapolation(
+ const bool & min,
+ const bool & max,
+ const IRtRelation & in_rt) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ RtParabolicExtrapolation rt_extrapolator; // r-t extrapolator
+ RtRelationLookUp *rt_low(0), *rt_high(0); // pointers to the r-t
+ // relationships after
+ // extrapolation
+ 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 = new 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 = new 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_low = new RtRelationLookUp(
+ rt_extrapolator.getRtWithParabolicExtrapolation(*rt_high,
+ 1.0,
+ 3.0,
+ 0.0,
+ add_fit_point));
+ } else {
+ rt_low = new RtRelationLookUp(
+ rt_extrapolator.getRtWithParabolicExtrapolation(in_rt,
+ 1.0,
+ 3.0,
+ 0.0,
+ add_fit_point));
+ }
+ }
+
+////////////////////
+// RETURN RESULTS //
+////////////////////
+
+ if (min && max) {
+ if(rt_high) delete rt_high;
+ if(in_rt.HasTmaxDiff())
+ {
+ rt_low->SetTmaxDiff(in_rt.GetTmaxDiff());
+ }
+ return rt_low;
+ }
+ if (min) {
+ if(rt_high) delete rt_high;
+ if(in_rt.HasTmaxDiff())
+ {
+ rt_low->SetTmaxDiff(in_rt.GetTmaxDiff());
+ }
+ return rt_low;
+ }
+ if(rt_low) delete rt_low;
+ if(in_rt.HasTmaxDiff())
+ {
+ rt_high->SetTmaxDiff(in_rt.GetTmaxDiff());
+ }
+ return rt_high;
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationIntegration.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationIntegration.cxx
new file mode 100644
index 0000000000000..7955289a61546
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtCalibrationIntegration.cxx
@@ -0,0 +1,416 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 02.02.2007, AUTHOR: OLIVER KORTNER
+// 12.11.2008, JOERG v. LOEBEN:
+// line 100: add missing if statement for close hits
+// line 140-240: change bining and adjust start/end points
+// line 252: change to r-t chebychev. this method does interpolating
+// the samplepoints better. the autocalibration methods are
+// working more accurate.
+// line 256: implement parabolic extrapolation for large driftradii.
+// 11.02.2010, include r(t_max)=r_max as addtional fit-point in extrapolation
+//
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS RtCalibrationIntegration ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// standard C++ //
+#include <algorithm>
+#include <iostream>
+#include <fstream>
+
+// MuonCalib //
+#include "MdtCalibRt/RtCalibrationIntegration.h"
+#include "MdtCalibT0/T0MTHistos.h"
+#include "MdtCalibData/RtFromPoints.h"
+#include "MdtCalibRt/RtParabolicExtrapolation.h"
+
+//root//
+#include "TF1.h"
+#include "TList.h"
+//::::::::::::::::::::::::
+//:: NAMESPACE SETTINGS ::
+//::::::::::::::::::::::::
+
+using namespace std;
+using namespace MuonCalib;
+
+inline Double_t slope_function_C(Double_t *x, Double_t *par)
+ {
+ Double_t &t(x[0]);
+ Double_t &a(par[0]),
+ &b(par[1]),
+ &t_0(par[2]),
+ &back(par[3]);
+ return back + exp(a+b*(t-t_0));
+ }
+
+inline void update_parameter_on_mttmax(TH1 * h, TF1 *f, const float &b, const float & T, const T0MTSettingsTMax &tmax_settings)
+ {
+ Double_t rmin, rmax;
+ f->GetRange(rmin, rmax);
+ TF1 *slope_function=new TF1("slope_function", slope_function_C, rmin, rmin+tmax_settings.WidthAB(), 4);
+ slope_function->SetParameter(0, f->GetParameter(T0MTHistos::TMAX_PAR_NR_A));
+ slope_function->FixParameter(1, b);
+ slope_function->FixParameter(2, f->GetParameter(T0MTHistos::TMAX_PAR_NR_T0));
+ slope_function->FixParameter(3, f->GetParameter(T0MTHistos::TMAX_PAR_NR_BACK));
+ h->Fit("slope_function", "R+", "");
+ f->FixParameter(T0MTHistos::TMAX_PAR_NR_A, slope_function->GetParameter(0));
+ f->FixParameter(T0MTHistos::TMAX_PAR_NR_B, b);
+ f->FixParameter(T0MTHistos::TMAX_PAR_NR_T, T);
+// std::cout<<"LLllLL "<<slope_function->GetParameter(0)<<" "<<b<<std::endl;
+// std::cout<<"LLllLL "<<f->GetParameter(T0MTHistos::TMAX_PAR_NR_A)<<" "<<f->GetParameter(T0MTHistos::TMAX_PAR_NR_B)<<std::endl;
+
+ }
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD init ::
+//:::::::::::::::::
+
+void RtCalibrationIntegration::init(bool close_hits,
+ double r_max,
+ double lower_extrapolation_radius,
+ double upper_extrapolation_radius,
+ bool add_tmax_difference){
+
+ m_close_hits = close_hits;
+ m_rt = 0;
+ m_r_max = r_max;
+ m_lower_extrapolation_radius = lower_extrapolation_radius;
+ m_upper_extrapolation_radius = upper_extrapolation_radius;
+ m_output = 0;
+ m_nb_hits_used = 0;
+ m_nb_segments_used = 0;
+ m_add_tmax_difference = add_tmax_difference;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::
+//:: METHOD number_of_hits_used ::
+//::::::::::::::::::::::::::::::::
+
+unsigned int RtCalibrationIntegration::number_of_hits_used(void) const {
+
+ return m_nb_hits_used;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::
+//:: METHOD analyseSegments ::
+//::::::::::::::::::::::::::::
+
+const IMdtCalibrationOutput
+*RtCalibrationIntegration::analyseSegments(const std::vector<MuonCalibSegment*> & seg) {
+
+ for (unsigned int k=0; k<seg.size(); k++) {
+ handleSegment(*seg[k]);
+ }
+ m_nb_segments_used = seg.size();
+ analyse();
+
+ return getResults();
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD handleSegment ::
+//::::::::::::::::::::::::::
+
+bool RtCalibrationIntegration::handleSegment(MuonCalibSegment & seg) {
+
+ ///////////////////////////////////////////////////////
+ // LOOP OVER THE HITS OF THE SEGMENTS AND STORE THEM //
+ ///////////////////////////////////////////////////////
+
+ // start with hits on the segment //
+ for (unsigned int k=0; k<seg.mdtHitsOnTrack(); k++) {
+ if(seg.mdtHOT()[k]->driftTime() <-8e8) continue;
+ m_t_drift.push_back(std::pair<double, bool>(seg.mdtHOT()[k]->driftTime(), seg.mdtHOT()[k]->identify().mdtMultilayer()==2));
+ }
+
+
+ // continue with close hits if requested //
+ if (m_close_hits==true){
+ for (unsigned int k=0; k<seg.mdtCloseHits(); k++) {
+ m_t_drift.push_back(std::pair<double, bool>(seg.mdtHOT()[k]->driftTime(), seg.mdtHOT()[k]->identify().mdtMultilayer()==2));
+ }
+ }
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD setInput ::
+//:::::::::::::::::::::
+
+void
+RtCalibrationIntegration::setInput(const IMdtCalibrationOutput * /*rt_input*/) {
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::
+//:: METHOD analyse ::
+//::::::::::::::::::::
+
+bool RtCalibrationIntegration::analyse(void) {
+
+ ///////////////
+ // VARIABLES //
+ ///////////////
+
+ T0MTSettings t0_setting; // settings of the MT t0 fitter
+ t0_setting.AddFitfun()=true;
+// t0_setting.DrawDebugGraphs()=true;
+ t0_setting.T0Settings()->ScrambleThreshold()=2;
+ t0_setting.T0Settings()->SlicingThreshold()=3;
+ t0_setting.TMaxSettings()->DistAB()+=50;
+ T0MTHistos drift_time_spec; // drift time spectrum used in the t0 and
+ // the tmax fit
+ T0MTHistos drift_time_spec_ml[2];
+
+
+ double t0, tmax; // t0 and tmax
+ int k_min(-1);//, k_max(-1); // first and last drift-time entry to be
+ // used in the integration procedure
+ unsigned int nb_bins(100); // number of integration bins
+ double bin_content; // number of entries per bin
+ vector<SamplePoint> point(nb_bins+1); // (t, r) points
+ double radius(0.0); // r(t)
+ double scf; // scale factor (r_max/number of used hits)
+ RtFromPoints rt_from_points; // r-t converter
+
+ ///////////////////////////////////////////////////////////////////////////
+ // STOP HERE, IF THERE ARE NOT ENOUGH ENTRIES IN THE DRIFT-TIME SPECTRUM //
+ ///////////////////////////////////////////////////////////////////////////
+
+ if (m_t_drift.size()<2000) {
+ cerr << endl
+ << "Class RtCalibrationIntegration, method analyse: "
+ << "FAILURE!\n"
+ << "Less than 2000 drift-time entries! No r-t "
+ << "relationship will be determined!\n";
+ return false;
+ }
+
+ ////////////////////////
+ // INTEGRATION METHOD //
+ ////////////////////////
+
+ // sort the hits in ascending order of their drift times //
+ sort(m_t_drift.begin(), m_t_drift.end());
+
+ // perform a t0 and tmax fit //
+ // fill the histogram with the drift-time spectrum //
+ int n_bins=static_cast<int>(32 * (200.0 + m_t_drift[m_t_drift.size() - 1].first - m_t_drift[0].first) / 25.0);
+ float min_t=m_t_drift[0].first - 100.0;
+ float max_t=m_t_drift[m_t_drift.size() - 1].first + 100.0;
+ TH1F *tspec= new TH1F("tspec", "DRIFT-TIME SPECTRUM", n_bins, min_t, max_t);
+ TH1F *tspec_ml[2];
+ tspec_ml[0]=new TH1F("tspec_ml0", "DRIFT-TIME SPECTRUM ML 0", n_bins, min_t, max_t);
+ tspec_ml[1]=new TH1F("tspec_ml1", "DRIFT-TIME SPECTRUM ML 1", n_bins, min_t, max_t);
+
+
+ for (unsigned int k=0; k<m_t_drift.size(); k++) {
+ tspec->Fill(m_t_drift[k].first, 1.0);
+ tspec_ml[static_cast<unsigned int>(m_t_drift[k].second)]->Fill(m_t_drift[k].first, 1.0);
+ }
+ drift_time_spec.SetTSpec(1, tspec, &t0_setting, false);
+ drift_time_spec_ml[0].SetTSpec(2, tspec_ml[0], &t0_setting, false);
+ drift_time_spec_ml[1].SetTSpec(3, tspec_ml[1], &t0_setting, false);
+
+ // t0 and tmax fits //
+
+ if (!drift_time_spec.FitT0() || !drift_time_spec.T0Ok()) {
+ cerr << endl
+ << "Class RtCalibrationIntegration, method analyse: "
+ << "FAILURE!\n"
+ << "t0 fit not successful, no r-t relationship will be "
+ << "calculated!\n";
+ return false;
+ }
+ t0 = drift_time_spec.GetT0Function()->GetParameter(
+ T0MTHistos::T0_PAR_NR_T0);
+
+ if (!drift_time_spec.FitTmax() || !drift_time_spec.TmaxOk()) {
+ cerr << endl
+ << "Class RtCalibrationIntegration, method analyse: "
+ << "FAILURE!\n"
+ << "tmax fit not successful, no r-t relationship will "
+ << "be calculated!\n";
+ return false;
+ }
+ tmax = drift_time_spec.GetTMaxFunction()->GetParameter(
+ T0MTHistos::TMAX_PAR_NR_TMAX)
+ +2.0*drift_time_spec.GetTMaxFunction()->GetParameter(
+ T0MTHistos::TMAX_PAR_NR_T);
+
+ // determine (t,r) points by integration //
+ m_nb_hits_used = 0;
+ for (unsigned int k=0; k<m_t_drift.size(); k++) {
+ if (m_t_drift[k].first>=t0 && m_t_drift[k].first<=tmax) {
+ m_nb_hits_used++;
+ }
+ if (k_min<0 && m_t_drift[k].first>=t0) {
+ k_min = k;
+ }
+ }
+ //k_max = k_min+m_nb_hits_used-1;
+
+ bin_content = static_cast<double>(m_nb_hits_used)/static_cast<double>(nb_bins);
+
+ scf = m_r_max/static_cast<double>(m_nb_hits_used);
+
+ point[0].set_x1(t0);
+ point[0].set_x2(0.0);
+ point[0].set_error(0.1); // give a higher weight to the end point in
+
+ // the final fit
+ for (unsigned int k=1; k<nb_bins; k++) {
+ radius = radius+scf*bin_content;
+ point[k].set_x1(m_t_drift[k_min+static_cast<int>(bin_content)*(k)].first );
+ point[k].set_x2(radius);
+ point[k].set_error(1.0);
+ }
+
+ point[nb_bins].set_x1(tmax);
+ point[nb_bins].set_x2(m_r_max);
+ point[nb_bins].set_error(1.);
+
+
+ /*//dump points
+ for(unsigned int i=0; i<point.size(); i++)
+ {
+ std::cout<<i<<" "<<point[i].x1()<<" "<<point[i].x2()<<" "<<point[i].error()<<std::endl;
+ }*/
+
+ // get the r-t relationship //
+ m_rt = new RtChebyshev(rt_from_points.getRtChebyshev(point, 15));
+ // m_rt = new RtRelationLookUp(rt_from_points.getRtRelationLookUp(point));
+
+ ///////////////////////////////////////////////////
+ // PARABOLIC EXTRAPOLATION FOR LARGE DRIFT RADII //
+ ///////////////////////////////////////////////////
+ vector<SamplePoint> add_fit_point; // additional r-t points for r(t_max)
+ add_fit_point.push_back(SamplePoint(tmax,m_r_max, 1.0));
+ RtParabolicExtrapolation rt_extrapolated;
+ RtRelationLookUp tmp_rt( rt_extrapolated.getRtWithParabolicExtrapolation(
+ *m_rt,
+ m_lower_extrapolation_radius,
+ m_upper_extrapolation_radius,
+ m_r_max,
+ add_fit_point));
+ IRtRelation *rt_new(0);
+ rt_new = new RtRelationLookUp(tmp_rt);
+
+
+ //////////////////////////////////////////////////
+ // Get length difference between multilayers //
+ //////////////////////////////////////////////////
+
+ if(tspec_ml[0]->GetEntries()>=10000 && tspec_ml[0]->GetEntries()>10000)
+ {
+ bool fit_ok(true);
+ float b[2];
+ float tmax[2];
+ float T[2];
+ for(unsigned int i=0; i<2; i++)
+ {
+ if(!drift_time_spec_ml[i].FitT0())
+ {
+ fit_ok=false;
+ break;
+ }
+ if(!drift_time_spec_ml[i].FitTmax())
+ {
+ fit_ok=false;
+ break;
+ }
+ b[i] = drift_time_spec_ml[i].GetTMaxFunction()->GetParameter(T0MTHistos::TMAX_PAR_NR_B);
+ tmax[i] = drift_time_spec_ml[i].GetTMaxFunction()->GetParameter(T0MTHistos::TMAX_PAR_NR_TMAX) - drift_time_spec_ml[i].GetT0Function()->GetParameter(T0MTHistos::T0_PAR_NR_T0);
+ T[i] = drift_time_spec_ml[i].GetTMaxFunction()->GetParameter(T0MTHistos::TMAX_PAR_NR_T);
+ }
+ if(fit_ok)
+ {
+ int refit=static_cast<int>((b[1] + 1.33717e-03)>(b[0] + 1.33717e-03));
+ int norefit=static_cast<bool>(refit)?0:1;
+// std::cout<<"HHhhHH "<<refit<<" "<<norefit<<std::endl;
+// std::cout<<"HHhhHH "<<b[norefit]<<std::endl;
+ TF1 *fixfun=drift_time_spec_ml[refit].GetTMaxFunctionNC();
+ fixfun->FixParameter(T0MTHistos::TMAX_PAR_NR_B, b[norefit]);
+ update_parameter_on_mttmax(drift_time_spec_ml[refit].GetTSpec(), fixfun, b[norefit], T[norefit],*t0_setting.TMaxSettings());
+ TList *l=drift_time_spec_ml[refit].GetTSpec()->GetListOfFunctions();
+ l->Remove(l->FindObject("mt_tmax_fermi"));
+ fit_ok=drift_time_spec_ml[refit].FitTmax();
+ tmax[refit] = drift_time_spec_ml[refit].GetTMaxFunction()->GetParameter(T0MTHistos::TMAX_PAR_NR_TMAX) - drift_time_spec_ml[refit].GetT0Function()->GetParameter(T0MTHistos::T0_PAR_NR_T0);
+ }
+
+ if(fit_ok && m_add_tmax_difference)
+ rt_new->SetTmaxDiff(tmax[0] - tmax[1]);
+ }
+
+
+ //////////////////////////////////////////////////
+ // STORE THE RESULTS IN THE RtCalibrationOutput //
+ //////////////////////////////////////////////////
+
+ if (m_output!=0) {
+ delete m_output;
+ }
+ m_output = new RtCalibrationOutput(rt_new,
+ new RtFullInfo("RtCalibrationIntegration", 1,
+ m_nb_segments_used, 0.0, 0.0, 0.0, 0.0));
+
+
+ return true;
+
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD converged ::
+//::::::::::::::::::::::
+
+bool RtCalibrationIntegration::converged(void) const {
+
+ return (m_output!=0);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD getResults ::
+//:::::::::::::::::::::::
+
+const IMdtCalibrationOutput * RtCalibrationIntegration::getResults(void) const {
+
+ return static_cast<const IMdtCalibrationOutput *>(m_output);
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtParabolicExtrapolation.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtParabolicExtrapolation.cxx
new file mode 100644
index 0000000000000..90ab4b98359f3
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibRt/src/RtParabolicExtrapolation.cxx
@@ -0,0 +1,248 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 10.11.2008, AUTHOR: OLIVER KORTNER
+// Modified: 28.02.2009 by O. Kortner and J. von Loeben, extend extrapolation
+// features
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MuonCalibMath/BaseFunctionFitter.h"
+#include "MuonCalibMath/LegendrePolynomial.h"
+#include "MdtCalibRt/RtParabolicExtrapolation.h"
+#include "MdtCalibData/IRtRelation.h"
+#include "MdtCalibData/RtRelationLookUp.h"
+
+//::::::::::::::::::::::::
+//:: NAMESPACE SETTINGS ::
+//::::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS RtCalibrationAnalytic ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: DEFAULT CONSTRUCTOR ::
+//:::::::::::::::::::::::::
+
+RtParabolicExtrapolation::RtParabolicExtrapolation(void) {
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD getRtWithParabolicExtrapolation(.,.,.) ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::
+
+RtRelationLookUp RtParabolicExtrapolation::getRtWithParabolicExtrapolation(
+ const IRtRelation & in_rt,
+ const double & r_min,
+ const double & r_max) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ double t_min(t_from_r(r_min, in_rt));
+ double t_max(t_from_r(r_max, in_rt));
+ vector<SamplePoint> t_r(10);
+ BaseFunctionFitter fitter(3);
+ LegendrePolynomial pol;
+ vector<double> rt_param(102); // parameters of the output r-t relationship
+ double r, t; // drift radius
+
+////////////////////////
+// FILL SAMPLE POINTS //
+////////////////////////
+
+ double step_size((t_max-t_min)/static_cast<double>(t_r.size()-1));
+ for (unsigned int k=0; k<t_r.size(); k++) {
+ t_r[k].set_x1(t_min+k*step_size);
+ t_r[k].set_x2(in_rt.radius(t_min+k*step_size));
+ t_r[k].set_error(1.0);
+ }
+ fitter.fit_parameters(t_r, 1, t_r.size(), &pol);
+
+ rt_param[0] = in_rt.tLower();
+ rt_param[1] = (in_rt.tUpper()-in_rt.tLower())/99.0;
+ for (unsigned int k=0; k<100; k++) {
+ t = rt_param[0]+rt_param[1]*k;
+ r = in_rt.radius(t);
+ if (r<r_min) {
+ rt_param[k+2] = r;
+ } else {
+ rt_param[k+2] = 0.0;
+ for (unsigned int l=0; l<3; l++) {
+ rt_param[k+2] = rt_param[k+2]+
+ fitter.coefficients()[l]*pol.value(l, t);
+ }
+ }
+ }
+
+ return RtRelationLookUp(rt_param);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD getRtWithParabolicExtrapolation(.,.,.,.,.) ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+RtRelationLookUp RtParabolicExtrapolation::getRtWithParabolicExtrapolation(
+ const IRtRelation & in_rt,
+ const double & r_min,
+ const double & r_max,
+ const double & r_ext,
+ const vector<SamplePoint> & add_fit_points
+ ) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+// double t_min(t_from_r(r_min, in_rt));
+ double t_min(get_max_t_at_r(r_min, in_rt));
+ double t_max(t_from_r(r_max, in_rt));
+ vector<SamplePoint> t_r(10);
+ BaseFunctionFitter fitter(3);
+ LegendrePolynomial pol;
+ vector<double> rt_param(102); // parameters of the output r-t relationship
+ double r, t; // drift radius, drift time
+
+////////////////////////
+// FILL SAMPLE POINTS //
+////////////////////////
+
+ double step_size((t_max-t_min)/static_cast<double>(t_r.size()-1));
+
+// r-t-points in fit region
+ for (unsigned int k=0; k<t_r.size(); k++) {
+ t_r[k].set_x1(t_min+k*step_size);
+ t_r[k].set_x2(in_rt.radius(t_min+k*step_size));
+ t_r[k].set_error(1.0);
+ }
+
+// addtional points for the fit
+ for (unsigned int k=0; k<add_fit_points.size(); k++){
+ t_r.push_back(add_fit_points[k]);
+ }
+
+// perform fit //
+ fitter.fit_parameters(t_r, 1, t_r.size(), &pol);
+
+// bring r-t-points in the right format for RtRelationLookUp and fill non
+// modified and extrapolated points in the r-t
+ rt_param[0] = in_rt.tLower();
+ rt_param[1] = (in_rt.tUpper()-in_rt.tLower())/99.0;
+// cerr <<"r_ext:\t"<< r_ext << "\tr_min\t" << r_min << "\tr_max\t" << r_max << endl;
+ for (unsigned int k=0; k<100; k++) {
+ t = rt_param[0]+rt_param[1]*k;
+ r = in_rt.radius(t);
+
+// distinguish if extrapolation area is right or left from fit region
+ if (r_ext<r_min) {
+ if (r>r_min && t>t_min) {
+ rt_param[k+2] = r; //fill original values
+ } else {
+ rt_param[k+2] = 0.0;
+ for (unsigned int l=0; l<3; l++) { // fill extrapolated values
+ rt_param[k+2] = rt_param[k+2]+
+ fitter.coefficients()[l]*pol.value(l, t);
+ }
+ if (rt_param[k+2]<0.0) {
+ rt_param[k+2] = 0.0;
+ }
+ }
+ }
+ else if (r_ext>r_max) {
+ if (r<r_max) {
+ rt_param[k+2] = r;//fill original values
+ } else {
+ rt_param[k+2] = 0.0;
+ for (unsigned int l=0; l<3; l++) { // fill extrapolated values
+ rt_param[k+2] = rt_param[k+2]+
+ fitter.coefficients()[l]*pol.value(l, t);
+ }
+ }
+ }
+// fill only original values if the radius where we want to extrapolate to is
+// within the fit region
+ else if (r_ext <= r_max && r_ext >=r_min) {
+ rt_param[k+2] = r;
+ cerr << "Class RtParabolicExtrapolation,"
+ << "WARNING: \n\t Extrapolated radius withing fit region"
+ << " - Nothing to be done."
+ << endl;
+
+ }
+ }
+
+ return RtRelationLookUp(rt_param);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD t_from_r ::
+//:::::::::::::::::::::
+
+double RtParabolicExtrapolation::t_from_r(const double & r,
+ const IRtRelation & in_rt) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ double precision(0.010); // spatial precision of the inversion
+ double t_max(in_rt.tUpper()); // upper time search limit
+ double t_min(in_rt.tLower()); // lower time search limit
+
+/////////////////////////////////////////////
+// SEARCH FOR THE CORRESPONDING DRIFT TIME //
+/////////////////////////////////////////////
+
+ while (t_max-t_min>0.1 &&
+ fabs(in_rt.radius(0.5*(t_min+t_max))-r)>precision) {
+
+ if (in_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);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::
+//:: METHOD get_max_t_at_r ::
+//:::::::::::::::::::::::::::
+
+double RtParabolicExtrapolation::get_max_t_at_r(const double & r,
+ const IRtRelation & in_rt) const {
+
+ for (double t=in_rt.tUpper(); t>=in_rt.tLower(); t=t-1.0) {
+ if (in_rt.radius(t)<r) {
+ return t+1.0;
+ }
+ }
+
+ return in_rt.tLower();
+
+}
--
GitLab