diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedCandidateFinder.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedCandidateFinder.h
new file mode 100644
index 0000000000000000000000000000000000000000..150b8d574891986b20df89644ed373d969d25792
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedCandidateFinder.h
@@ -0,0 +1,73 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 13.07.2008, AUTHOR: OLIVER KORTNER
+// Modified: 04.08.2008 by O. Kortner, estimated direction of incidence can be
+// set.
+// 07.08.2008 by O. Kortner, bug fig in the pattern recognition.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef MuonCalib_CurvedCandidateFinderH
+#define MuonCalib_CurvedCandidateFinderH
+
+//:::::::::::::::::::::::::::::::::
+//:: CLASS CurvedCandidateFinder ::
+//:::::::::::::::::::::::::::::::::
+
+/// \class CurvedCandidateFinder
+///
+/// This class searches for candidates of curved parabolic lines connecting
+/// all hits given to the class within a user defined search road.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 13.07.2008
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// STL //
+#include <vector>
+
+// MuonCalib //
+#include "MuonCalibEventBase/MdtCalibHitBase.h"
+#include "MdtCalibFitters/CurvedLine.h"
+
+
+namespace MuonCalib {
+
+class CurvedCandidateFinder {
+
+public:
+// Constructors //
+ CurvedCandidateFinder(const std::vector<const MdtCalibHitBase *> & hits);
+ ///< Constructor
+ ///< @param hits Vector of hits used in the candidate finding.
+
+// Methods //
+ const std::vector<CurvedLine> & getCandidates(const double & road_width);
+ ///< get all candidates connecting all hits
+ ///< within the given road width (mm)
+ const std::vector<CurvedLine> & getCandidates(const double & road_width,
+ const Amg::Vector3D & est_dir);
+ ///< get all candidates connecting all hits
+ ///< within the given road width (mm);
+ ///< est_dir is the estimated direction of
+ ///< incidence
+
+private:
+// hits //
+ std::vector<const MdtCalibHitBase *> m_hits; // vector hits used in the
+ // candidate search
+
+// candidates //
+ std::vector<CurvedLine> m_candidates; // vector of candidate lines
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedLine.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedLine.h
new file mode 100644
index 0000000000000000000000000000000000000000..9177eead1449407bdef97465d35ee16fb437d4c0
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedLine.h
@@ -0,0 +1,93 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 05.04.2008, AUTHOR: OLIVER KORTNER
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef MuonCalib_CurvedLineH
+#define MuonCalib_CurvedLineH
+
+//::::::::::::::::::::::
+//:: CLASS CurvedLine ::
+//::::::::::::::::::::::
+
+/// \class CurvedLine
+///
+/// This class describes a curved line in an MDT chamber. Let x,y,z denote the
+/// local coordinate frame of an MDT chamber. Then the class assumes a
+/// straight line in the xz and a parabolicly curved line in the yz plane.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 05.04.2008
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// STL //
+#include <vector>
+
+// CLHEP //
+#include "MdtCalibFitters/MTStraightLine.h"
+
+// MuonCalib //
+#include "MuonCalibMath/Legendre_polynomial.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "EventPrimitives/EventPrimitives.h"
+
+namespace MuonCalib {
+
+class CurvedLine {
+
+public:
+// Constructors //
+ CurvedLine(void);
+ ///< Default constructor: a straight line through (0,0,0) pointing in
+ ///< in the local x direction of the chambers is created.
+
+ CurvedLine(std::vector<Amg::Vector3D> & points);
+ ///< Constructor.
+ ///< \param points A straight line in the xz and a parabolic line in the
+ ///< yz plane is fitter through the given points.
+ ///< All points get equal weight.
+
+ CurvedLine(std::vector<Amg::Vector3D> & points,
+ std::vector<Amg::Vector3D> & x_and_y_errors);
+ ///< Constructor.
+ ///< \param points A straight line in the xz and a parabolic line in the
+ ///< yz plane is fitter through the given points.
+ ///< \param x_and_y_errors Errors of the x and y coordinates of the given
+ ///< points
+
+// Methods //
+// get-methods //
+ Amg::Vector3D getPointOnLine(const double & loc_z) const;
+ ///< get the point on the line a the local
+ ///< z coordinate "loc_z"
+
+ MTStraightLine getTangent(const double & loc_z) const;
+ ///< get the tangent to the line a the local
+ ///< z coordinate "loc_z"
+
+private:
+// parameters of the curved line //
+ Legendre_polynomial *m_Legendre; // pointer to the Legendre polynomial
+ // describing the curved line
+ Amg::VectorX m_coeff_xz; // coefficients of the straight line in the local
+ // xz plane
+ Amg::VectorX m_coeff_yz; // coefficients of the curved line in the local
+ // yz plane
+
+// private methods //
+ void init(std::vector<Amg::Vector3D> & points,
+ std::vector<Amg::Vector3D> & x_and_y_errors);
+ // initialization routine
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedPatRec.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedPatRec.h
new file mode 100644
index 0000000000000000000000000000000000000000..3f0b056d6db78445d6cbbe8e68270025b14444b3
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/CurvedPatRec.h
@@ -0,0 +1,172 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 05.04.2008, AUTHOR: OLIVER KORTNER
+// Modified: 25.07.2008 by O. Kortner, improved pattern recognition by using
+// the class "CurvedCandidateFinder".
+// 04.08.2008 by O. Kortner, further improvements of the pattern
+// recognition for large incidence angles.
+// 07.08.2008 by O. Kortner, bug fig in the pattern recognition.
+// 18.08.2008 by O. Kortner, update of chi^2 and segment position
+// and direction added.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef MuonCalib_CurvedPatRecH
+#define MuonCalib_CurvedPatRecH
+
+//::::::::::::::::::::::::
+//:: CLASS CurvedPatRec ::
+//::::::::::::::::::::::::
+
+/// \class CurvedPatRec
+///
+/// This class searches for the best hit combination for a parabolic curved
+/// segment fit and fits a parabola in the precision plane to the hits found
+/// before. The algorithm overwrites the segment position nor the
+/// segment direction as well as the chi^2 in the latest versions.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 18.08.2008
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "GeoPrimitives/GeoPrimitives.h"
+
+// STL //
+#include <vector>
+
+// MuonCalib //
+#include "MdtCalibInterfaces/IMdtPatRecFitter.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "MdtCalibFitters/MTStraightLine.h"
+#include "MdtCalibFitters/StraightPatRec.h"
+#include "MdtCalibFitters/CurvedLine.h"
+
+
+namespace MuonCalib {
+
+class CurvedPatRec : public IMdtPatRecFitter {
+
+public:
+// Constructors //
+ CurvedPatRec(void);
+ ///< Default constructor: road width of 0.5 mm is used.
+
+ CurvedPatRec(const double & road_width);
+ ///< Constructor: user defined road width [mm] for pattern recognition.
+
+// Methods //
+ double roadWidth(void) const;
+ ///< get the road width used in the
+ ///< pattern recognition [mm]
+ unsigned int numberOfTrackHits(void) const;
+ ///< get the number of track hits in the
+ ///< final track,
+ ///< the method returns 0, if no track
+ ///< has been reconstructed
+ const std::vector<const MdtCalibHitBase*> & trackHits(void) const;
+ ///< get a vector with a pointer to the
+ ///< hits used in the track
+ ///< reconstruction
+ double chi2(void) const;
+ ///< get the chi^2 of the final track,
+ ///< the method returns -1, if no track
+ ///< has been reconstructed
+ double chi2PerDegreesOfFreedom(void) const;
+ ///< get the chi^2 per degrees of
+ ///< freedom for the final track,
+ ///< the method returns -1, if no track
+ ///< has been reconstructed
+ const CurvedLine & curvedTrack(void) const;
+ ///< get the reconstructed curved track
+
+// set-method //
+ void setRoadWidth(const double & r_road_width);
+ ///< set the road width [mm] for the pattern
+ ///< recognition = r_road_width
+ void setTimeOut(const double & time_out);
+ ///< set the time-out for the track
+ ///< finding to time_out (in seconds)
+
+// methods required by the base class "IMdtSegmentFitter" //
+ bool fit(MuonCalibSegment & r_segment) const;
+ ///< reconstruction of the track using
+ ///< all hits in the segment
+ ///< "r_segment", returns true in case
+ ///< of success;
+ ///< the algorithm overwrites
+ ///< the track radii, the track
+ ///< position, track direction, and
+ ///< chi^2 per degrees of freedom;
+ ///< warning: the errors of the track
+ ///< radii are only approximate
+ bool fit(MuonCalibSegment & r_segment, HitSelection r_selection) const;
+ ///< reconstruction of the track using
+ ///< only those hits in r_segment for
+ ///< which the r_selection[.] is 0,
+ ///< return true in case of success;
+ ///< the algorithm overwrites
+ ///< the track position, direction,
+ ///< and the chi^2 in r_segment; it
+ ///< updates the distances of all hits
+ ///< from the track, i.e. also of those
+ ///< hits which were rejected from the
+ ///< track reconstruction;
+ ///< warning: the errors of the track
+ ///< radii are only approximate
+ void printLevel(int level);
+ ///< set the print level, this call
+ ///< has no effect
+
+private:
+// internal co-ordinate definition //
+// x3
+// ^
+// o o o o o o |
+// ... o o o o o ... o--> x2
+// o o o o o o x1
+//
+
+// settings //
+ double m_road_width; // road width to be used in pattern finding
+ double m_time_out; // time out for pattern recognition.
+
+// straight pattern recognition to get the approximate incidence angle //
+ mutable StraightPatRec m_sfitter;
+
+// final track parameters //
+ mutable std::vector<const MdtCalibHitBase*> m_track_hits;
+ // pointer to hits used by the final
+ // track
+ mutable double m_chi2; // chi^2 of the final track
+ mutable CurvedLine m_curved_track; // curved trajectory
+
+// auxiliary methods //
+ Amg::Vector3D getHitPoint(const MdtCalibHitBase * hit,
+ const MTStraightLine & straight_track) const;
+ // transform the given hit into a hit point in
+ // 3 dimensions using the given straight line
+ // as helper
+ std::vector<Amg::Vector3D> getHitPoints(
+ std::vector<const MdtCalibHitBase*> track_hits,
+ const MTStraightLine & straight_track) const;
+ // transform the track hits into hit points in
+ // 3 dimensions using the given straight line
+ // as helper
+ std::vector<Amg::Vector3D> getHitPoints(
+ std::vector<const MdtCalibHitBase*> track_hits,
+ const CurvedLine & curved_track) const;
+ // transform the track hits into hit points in
+ // 3 dimensions using the given curved line
+ // as helper
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h
new file mode 100644
index 0000000000000000000000000000000000000000..fc38a64e70c786a6fbfe800507e7da62c3950fbf
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/DCSLFitter.h
@@ -0,0 +1,63 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// DCSLFitter.h
+// Header file for class DCSLFitter
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+// nveldik@nikhef.nl
+///////////////////////////////////////////////////////////////////
+
+#ifndef MUONCALIB_DCSLFITTER_H
+#define MUONCALIB_DCSLFITTER_H
+
+#include "MdtCalibInterfaces/IMdtSegmentFitter.h"
+
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+
+namespace MuonCalib {
+/**@class DCSLFitter
+ Straight line fitter for drift circles
+
+ @author Niels.Van.Eldik@cern.ch
+*/
+
+
+ class DCSLFitter : public IMdtSegmentFitter {
+ public:
+ DCSLFitter() : m_debug(false), m_error_dy0(0.0), m_error_dtheta(0.0) {}
+
+ /** fit using all hits */
+ bool fit( MuonCalibSegment& seg ) const;
+
+ /** fit subset of the hits.
+ If the HitSelection vector contains a 0 for a given hit the hit is used
+ else the hit is not included in the fit.
+ The size of the HitSelection vector should be equal to the number of hits on track else
+ no fit is performed.
+ */
+ bool fit( MuonCalibSegment& seg, HitSelection selection ) const;
+// double error_dy0() { return m_error_dy0;}
+// double error_dtheta() { return m_error_dtheta;}
+
+ /** set print level */
+ void printLevel(int level);
+ private:
+
+ bool m_debug;
+
+ mutable double m_error_dy0;
+ mutable double m_error_dtheta;
+
+ double getY( const Amg::Vector3D& p ) const { return p.z(); }
+ double getZ( const Amg::Vector3D& p ) const { return p.y(); }
+ Amg::Vector3D getVec( double x, double y, double z) const { return Amg::Vector3D( x, z, y ); }
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/IndexSet.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/IndexSet.h
new file mode 100644
index 0000000000000000000000000000000000000000..4c4bf6f7aed19ae776a05068355f69d8ae57eb26
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/IndexSet.h
@@ -0,0 +1,104 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 09.05.2005, AUTHOR: OLIVER KORTNER
+// Modified: 16.07.2006 by O. Kortner, vector size check in initialization
+// added,
+// != operator corrected.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef IndexSetHXX
+#define IndexSetHXX
+
+//::::::::::::::::::::
+//:: CLASS IndexSet ::
+//::::::::::::::::::::
+
+/// \class IndexSet
+/// This class allows the user to save a set of indices and to compare two sets
+/// of indices. IndexSet is an auxiliary class for
+/// QuasianalyticLineReconstruction.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 09.05.2005
+
+//////////////////
+// HEADER FILES //
+//////////////////
+
+// standard C++ libraries //
+#include <iostream>
+
+// STEL //
+#include <vector>
+
+namespace MuonCalib {
+
+class IndexSet {
+
+private:
+// internal index memory //
+ unsigned int m_nb_indices; // number of indices
+ std::vector<int> m_index; // indices
+
+// initialization methods //
+ void m_init(void);
+ // default initiailization method; the number of indices is set to 0
+ void m_init(const unsigned int & r_nb_indices);
+ // initialization method; the number of indices is set to r_nb_indices
+ void m_init(const unsigned int & r_nb_indices,
+ const std::vector<int> r_index);
+ // initialization method; the number of indices is set to
+ // r_nb_indices, the vector r_index contains the indices
+
+public:
+// Constructors //
+ IndexSet(void) {
+ m_init();
+ }
+ ///< default constructor: the number of indices is set to 0
+
+ IndexSet(const unsigned int & r_nb_indices) {
+ m_init(r_nb_indices);
+ }
+ ///< constructor: the number of indices is set to r_nb_indices
+
+ IndexSet(const unsigned int & r_nb_indices,
+ const std::vector<int> r_index) {
+ m_init(r_nb_indices, r_index);
+ }
+ ///< constructor: the number of indices is set to r_nb_indices,
+ ///< the vector r_index contains the indices
+
+// Methods //
+ unsigned int size(void) const;
+ ///< get the number of indices
+ void resize(const unsigned int & r_size);
+ ///< resize the index set to r_size;
+ ///< the index store is preserved as
+ ///< far as possible
+ int operator [] (const unsigned int & r_k);
+ int operator [] (const unsigned int & r_k) const;
+ ///< access to the indices stored in
+ ///< the class,
+ ///< 0 <= r_k < size();
+ ///< WARNING: no test on the index is
+ ///< performed
+ void sort(void);
+ ///< sort the indices in ascending order
+ bool operator == (const IndexSet & r_index_set) const;
+ bool operator != (const IndexSet & r_index_set) const;
+ ///< comparison of two index sets;
+ ///< two index sets are considered
+ ///< equal if they are of the same size
+ ///< and their indices are the same at
+ ///< the same position [.]
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/LocalSegmentResolver.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/LocalSegmentResolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..cff59fbf4a584f91f583d2353cd45e2a4b58d7b4
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/LocalSegmentResolver.h
@@ -0,0 +1,59 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// LocalSegmentResolver.h
+// Header file for class MdtSegment
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+// nveldik@nikhef.nl
+///////////////////////////////////////////////////////////////////
+
+#ifndef MUONCALIB_LOCALSEGMENTRESOLVER_H
+#define MUONCALIB_LOCALSEGMENTRESOLVER_H
+
+
+// stl includes
+#include <vector>
+
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+
+class MdtCalibHitBase;
+
+namespace MuonCalib {
+
+/**@class LocalSegmentResolver
+ Class calculates the local position and direction of the segment using the hits on track.
+
+ @author Niels.Van.Eldik@cern.ch
+*/
+
+ class LocalSegmentResolver {
+ public:
+ /** constructor */
+ LocalSegmentResolver();
+
+ /** resolve local position and direction of the track segment */
+ bool resolve(MuonCalibSegment* seg) const;
+
+ /** set print level */
+ void setPrintLevel(int level) { m_printLevel = level; }
+
+ private:
+ typedef MuonCalibSegment::MdtHitVec HitVec;
+ typedef std::pair<Amg::Vector3D,Amg::Vector3D > Line;
+ typedef std::vector<Line> LineVec;
+ private:
+ LineVec getLines( const MdtCalibHitBase& firstHit, const MdtCalibHitBase& lastHit ) const;
+
+ int bestLine( const HitVec& hits, const LineVec& localTracks ) const;
+
+ /** print level */
+ int m_printLevel;
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/LocalToPrecision.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/LocalToPrecision.h
new file mode 100644
index 0000000000000000000000000000000000000000..f64660c461bc7ff9e19bffcf24682ade10e31100
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/LocalToPrecision.h
@@ -0,0 +1,22 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MUONCALIB_LOCALTOPRECISION_H
+#define MUONCALIB_LOCALTOPRECISION_H
+
+#include "GeoPrimitives/GeoPrimitives.h"
+
+namespace MuonCalib {
+
+ /** Provides transformation between local and presicion reference frames*/
+ class LocalToPrecision {
+ public:
+ static Amg::Vector3D precision( const Amg::Vector3D& p ) { return Amg::Vector3D( p.y(), p.z(), p.x() ); }
+ static Amg::Vector3D local( const Amg::Vector3D& p ) { return Amg::Vector3D( p.z(), p.x(), p.y() ); }
+
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/MTStraightLine.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/MTStraightLine.h
new file mode 100644
index 0000000000000000000000000000000000000000..ec2b3dc00518e8046183ecf3666d634e4cac05dc
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/MTStraightLine.h
@@ -0,0 +1,168 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 01.03.2006, AUTHORS: OLIVER KORTNER, FELIX RAUSCHER
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::::::::::
+//:: CLASS MtStraightLine ::
+//::::::::::::::::::::::::::
+#include "GeoPrimitives/GeoPrimitives.h"
+
+namespace MuonCalib {
+/// \class MtStraightLine
+/// This class contains a simple implementation of a straifht line with errors
+/// on the slope and the intercept of the straight line.
+///
+/// \author Oliver.Kortner@CERN.CH \author Felix.Rauscher@CERN.CH
+///
+/// \date 01.03.2006
+}
+
+#ifndef MTStraightLineH
+#define MTStraightLineH
+
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// CLHEP //
+//#include "CLHEP/config/CLHEP.h"
+//#include "CLHEP/Units/SystemOfUnits.h"
+//#include "CLHEP/Units/PhysicalConstants.h"
+
+namespace MuonCalib {
+
+class MTStraightLine {
+
+public:
+// Constructors //
+ MTStraightLine(void) {
+ m_init();
+ }
+ ///< Default constructor: all internal parameters of the straight line
+ ///< are set to 0.
+
+ MTStraightLine(const Amg::Vector3D & r_position,
+ const Amg::Vector3D & r_direction,
+ const Amg::Vector3D & r_position_error,
+ const Amg::Vector3D & r_direction_error) {
+ m_init(r_position, r_direction,
+ r_position_error, r_direction_error);
+ }
+ ///< Constructor:
+ ///< r_position: position vector of the straight line;
+ ///< r_direction: direction vector of the straight line;
+ ///< r_position_error: error on the position vector;
+ ///< r_direction_error: error on the direction vector.
+
+ MTStraightLine(const double & r_m_x1, const double & r_b_x1,
+ const double & r_m_x2, const double & r_b_x2,
+ const double & r_m_x1_err, const double & r_b_x1_err,
+ const double & r_m_x2_err, const double & r_b_x2_err) {
+ m_init(r_m_x1, r_b_x1, r_m_x2, r_b_x2,
+ r_m_x1_err, r_b_x1_err, r_m_x2_err, r_b_x2_err);
+ }
+ ///< Coordinates: x1, x2, x3.
+ ///< Parametrization of the straight line: x1 = r_m_x1*x3 + r_b_x1;
+ ///< x2 = r_m_x2*x3 + r_b_x2;
+ ///< x3 arbitrary.
+ ///< r_m_x1_err: error on r_m_x1.
+ ///< r_b_x1_err: error on r_b_x1.
+ ///< r_m_x2_err: error on r_m_x2.
+ ///< r_b_x2_err: error on r_b_x2.
+
+// Methods //
+// get-methods //
+ Amg::Vector3D positionVector(void) const;
+ ///< get the position vector of the
+ ///< straight line
+ Amg::Vector3D directionVector(void) const;
+ ///< get the direction vector of the
+ ///< straight line
+ Amg::Vector3D positionError(void) const;
+ ///< get the error on the position
+ ///< vector of the straight line
+ Amg::Vector3D directionError(void) const;
+ ///< get the error on the direction
+ ///< vector of the straight line
+ double m_x1(void) const;
+ ///< get the slope of the straight line
+ ///< in the x1-x3 plane
+ double m_x1_error(void) const;
+ ///< get the error on the slope of the
+ ///< straight line in the x1-x3 plane
+ double b_x1(void) const;
+ ///< get the intercept of the straight
+ ///< line in the x1-x3 plane
+ double b_x1_error(void) const;
+ ///< get the error on the intercept of
+ ///< the straight line in the x1-x3
+ ///< plane
+ double m_x2(void) const;
+ ///< get the slope of the straight line
+ ///< in the x2-x3 plane
+ double m_x2_error(void) const;
+ ///< get the error on the slope of the
+ ///< straight line in the x2-x3 plane
+ double b_x2(void) const;
+ ///< get the intercept of the straight
+ ///< line in the x2-x3 plane
+ double b_x2_error(void) const;
+ ///< get the slope of the intercept of
+ ///< the straight line in the x2-x3
+ ///< plane
+ Amg::Vector3D pointOnLine(const double & lambda) const;
+ ///< get the point on the line for the
+ ///< given scale factor lambda,
+ ///< point=position_vector+lambda*direction_vector
+ double signDistFrom(const MTStraightLine & h) const;
+ ///< get the signed distance of two
+ ///< lines
+ ///< (if both are parallel, dist>0)
+ double distFromLine(const Amg::Vector3D & point) const;
+ ///< get the distance of point point
+ ///< from straight line
+
+private:
+// internal representation of the straight line //
+ Amg::Vector3D m_position; //position vector of the straight line
+ Amg::Vector3D m_direction; //direction vector of the straight line
+
+// errors on the position and direction vector //
+ Amg::Vector3D m_position_error; //error on the position vector
+ Amg::Vector3D m_direction_error; //error on the direction vector
+
+// initializtion methods //
+ void m_init(void); //default initialization method
+
+ void m_init(const Amg::Vector3D & r_position,
+ const Amg::Vector3D & r_direction,
+ const Amg::Vector3D & r_position_error,
+ const Amg::Vector3D & r_direction_error);
+ // r_position: position vector of the straight line;
+ // r_direction: direction vector of the straight line;
+ // r_position_error: error on the position vector;
+ // r_direction_error: error on the direction vector
+
+ void m_init(const double & r_m_x1, const double & r_b_x1,
+ const double & r_m_x2, const double & r_b_x2,
+ const double & r_m_x1_err, const double & r_b_x1_err,
+ const double & r_m_x2_err, const double & r_b_x2_err);
+ // Coordinates: x1, x2, x3.
+ // Parametrization of the straight line: x1 = r_m_x1*x3 + r_b_x1;
+ // x2 = r_m_x2*x3 + r_b_x2;
+ // x3 arbitrary.
+ // r_m_x1_err: error on r_m_x1.
+ // r_b_x1_err: error on r_b_x1.
+ // r_m_x2_err: error on r_m_x2.
+ // r_b_x2_err: error on r_b_x2.
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/MuCCaFitter.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/MuCCaFitter.h
new file mode 100644
index 0000000000000000000000000000000000000000..67431848c8e988baf7a8f745608f8bab3d20426c
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/MuCCaFitter.h
@@ -0,0 +1,112 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+///////////////////////////////////////////////////////////////////
+// MuCCaFitter.h
+// Header file for class MuCCaFitter
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef MuCCaFITTER_H
+#define MuCCaFITTER_H
+
+#include "MdtCalibInterfaces/IMdtSegmentFitter.h"
+
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+
+#include <vector>
+
+
+namespace MuonCalib {
+
+/**@class MuCCaFitter
+ Interface to the straight line fitter for drift circles used by Calib
+
+ @author Fabrizio.Petrucci@cern.ch
+*/
+class MuCCaFitter : public IMdtSegmentFitter {
+ public:
+ /** fit using all hits */
+ bool fit( MuonCalibSegment& seg ) const;
+
+ /** fit subset of the hits.
+ If the HitSelection vector contains a 0 for a given hit the hit is used
+ else the hit is not included in the fit.
+ The size of the HitSelection vector should be equal to the number of hits on track else
+ no fit is performed.
+ */
+ bool fit( MuonCalibSegment& seg, HitSelection selection ) const;
+
+ /** set print level */
+ void printLevel(int level);
+
+ private:
+
+ bool m_debug;
+
+ /** these methods are needed to change the reference frame between the
+ local one of the hit and one used by the fitter: apparently the 2 frames
+ are the same ... TO BE CHECKED!!!
+ */
+ // double getY( const Amg::Vector3D& p ) const { return p.y(); }
+ // double getX( const Amg::Vector3D& p ) const { return p.x(); }
+ // double getZ( const Amg::Vector3D& p ) const { return p.z(); }
+ // Amg::Vector3D getVec( double x, double y, double z) const { return Amg::Vector3D( x, y, z ); }
+ double getY( const Amg::Vector3D& p ) const { return p.z(); }
+ double getX( const Amg::Vector3D& p ) const { return p.x(); }
+ double getZ( const Amg::Vector3D& p ) const { return p.y(); }
+ Amg::Vector3D getVec( double x, double y, double z) const { return Amg::Vector3D( x, z, y ); }
+};
+
+/**@class MuCCaFitterImplementation
+ The straight line fitter for drift circles used by Calib
+
+ @author Fabrizio.Petrucci@cern.ch, Paolo.Branchini@cern.ch
+*/
+class MuCCaFitterImplementation{
+ public:
+ void Computelin(double x1, double y1, double r1,
+ double x2, double y2, double r2);
+ void Computelinparnew(double x1, double y1, double r1,
+ double x2, double y2, double r2);
+ void computehitsfromcircles(double x0,double y0,double r0,
+ double x1,double y1,double r1,double a,double b);
+ // void Computehitpsemes(int nhit,double *xcirc, double *ycirc,
+ // double *rcirc, double a, double b);
+ void Computehitpsemes(int nhit,std::vector<double> xcirc, std::vector<double> ycirc, std::vector<double>rcirc, double a, double b);
+ // void Computeparam3(int number_of_hits,double x[100],double y[100],
+ // double r[100],double sr[100]);
+ void Computeparam3(int number_of_hits,std::vector<double> x,std::vector<double> y,std::vector<double> r,std::vector<double> sr);
+ double get_a();
+ double get_b();
+ double get_da();
+ double get_db();
+ double get_corrab();
+ double get_chi2f();
+ private:
+ /** parameters of the 4 tangent lines from first and last hit */
+ double m_angularcoefficient[4];
+ double m_bfpar[4];
+ /** track points from 2 circles */
+ double m_xout0;
+ double m_yout0;
+ double m_xout1;
+ double m_yout1;
+ /** track points */
+ double m_xpoint[100];
+ double m_ypoint[100];
+ /** results */
+ double m_aoutn; //!< track slope
+ double m_sig2a; //!< track slope's variance
+ double m_bout; //!< track constant term
+ double m_sig2b; //!< error on track constant term
+ double m_corrab;//!< correlation term
+ double m_chi2f; //!< chisquared
+};
+
+}
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/QuasianalyticLineReconstruction.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/QuasianalyticLineReconstruction.h
new file mode 100644
index 0000000000000000000000000000000000000000..02e1841ccca53cbf40b1b3576f8c8cdbf9ac6b15
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/QuasianalyticLineReconstruction.h
@@ -0,0 +1,233 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 01.03.2006, AUTHOR: OLIVER KORTNER
+// Modified: 15.07.2006 by O. Kortner, error calculation corrected,
+// chi^2 refit functionality added.
+// 13.01.2007 by O. Kortner, bug fix in candidate treatment, some
+// candidates were considered to be
+// identical although different;
+// modifications to improve the
+// reconstruction efficiency at very high
+// background count rates.
+// 27.03.2007 by O. Kortner, distances with signs filled into
+// MuonCalibSegment.
+// 23.03.2007 by O. Kortner, isnan check for chi^2.
+// 08.06.2007 by O. Kortner, final track segment has rphi position and
+// direction of the initial segment.
+// 23.06.2006 by O. Kortner, add convention for rphi track if the
+// pattern recognition has failed to
+// provide it.
+// 26.11.2007 by O. Kortner, fix for segment refinement.
+// 13.12.2007 by O. Kortner, time-out added.
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef QuasianalyticLineReconstructionH
+#define QuasianalyticLineReconstructionH
+
+//:::::::::::::::::::::::::::::::::::::::::::
+//:: CLASS QuasianalyticLineReconstruction ::
+//:::::::::::::::::::::::::::::::::::::::::::
+
+/// \class QuasianalyticLineReconstruction
+/// This class performs the reconstruction of straight line in MDT chambers by
+/// calculating tangents to all pairs of two tubes on a straight-line track and
+/// averaging over the calculated tangents to get the parameters of the
+/// straight-line track. The class redoes the pattern recognition using the
+/// hits provided to it. The hits which were used in the track fit can be
+/// retrieved from the class. A refit of the result of the quasianalytic
+/// reconstruction with the DCSLFitter can requested.
+/// The user can set a time-out for the track finding. It is set to 10 seconds
+/// by default.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 01.03.2006, 16.07.2006, 13.12.2007
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// CLHEP //
+//#include "CLHEP/config/CLHEP.h"
+
+// STL //
+#include <vector>
+
+// MuonCalib //
+#include "MdtCalibInterfaces/IMdtPatRecFitter.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "MdtCalibFitters/MTStraightLine.h"
+#include "MdtCalibFitters/IndexSet.h"
+#include "MdtCalibFitters/DCSLFitter.h"
+
+
+namespace MuonCalib {
+
+class QuasianalyticLineReconstruction : public IMdtPatRecFitter {
+
+public:
+// Constructors //
+ QuasianalyticLineReconstruction(void) {
+ m_init();
+ }
+ ///< Default constructor: road width for pattern recognition = 0.5 mm.
+
+ QuasianalyticLineReconstruction(const double & r_road_width) {
+ m_init(r_road_width);
+ }
+ ///< Constructor: user-defined road width for pattern recognition.
+
+// Methods //
+// get-methods //
+ double roadWidth(void) const;
+ ///< get the road width used in the
+ ///< pattern recognition
+ unsigned int numberOfTrackHits(void) const;
+ ///< get the number of track hits in the
+ ///< final track,
+ ///< the method returns 0, if no track
+ ///< has been reconstructed
+ const std::vector<const MdtCalibHitBase*> & trackHits(void) const;
+ ///< get a vector with a pointer to the
+ ///< hits used in the track
+ ///< reconstruction
+ double chi2(void) const;
+ ///< get the chi^2 of the final track,
+ ///< the method returns -1, if no track
+ ///< has been reconstructed
+ double chi2PerDegreesOfFreedom(void) const;
+ ///< get the chi^2 per degrees of
+ ///< freedom for the final track,
+ ///< the method returns -1, if no track
+ ///< has been reconstructed
+ MTStraightLine track(void) const;
+ ///< get the final track in the local
+ ///< co-ordinate frame, i.e.
+ ///< z
+ ///< ^
+ ///< o o o o o o |
+ ///< ... o o o o o ... o--> y
+ ///< o o o o o o x
+
+// set-method //
+ void setRoadWidth(const double & r_road_width);
+ ///< set the road width for the pattern
+ ///< recognition = r_road_width
+ void setTimeOut(const double & time_out);
+ ///< set the time-out for the track
+ ///< finding to time_out (in seconds)
+// methods required by the base class "IMdtSegmentFitter" //
+ bool fit(MuonCalibSegment & r_segment) const;
+ ///< reconstruction of the track using
+ ///< all hits in the segment
+ ///< "r_segment", returns true in case
+ ///< of success;
+ ///< the algorithm overwrites
+ ///< the track radii, the track
+ ///< position, track direction, and
+ ///< chi^2 per degrees of freedom;
+ ///< warning: the errors of the track
+ ///< radii are only approximate
+ bool fit(MuonCalibSegment & r_segment, HitSelection r_selection) const;
+ ///< reconstruction of the track using
+ ///< only those hits in r_segment for
+ ///< which the r_selection[.] is 0,
+ ///< return true in case of success;
+ ///< the algorithm overwrites
+ ///< the track position, direction,
+ ///< and the chi^2 in r_segment; it
+ ///< updates the distances of all hits
+ ///< from the track, i.e. also of those
+ ///< hits which were rejected from the
+ ///< track reconstruction;
+ ///< warning : the errors of the track
+ ///< radii are only approximate
+ void printLevel(int level);
+ ///< set the print level, this call
+ ///< has no effect
+
+private:
+// internal co-ordinate definition //
+// x3
+// ^
+// o o o o o o |
+// ... o o o o o ... o--> x2
+// o o o o o o x1
+//
+
+// final track parameters //
+ mutable std::vector<const MdtCalibHitBase*> m_track_hits;
+ // pointer to hits used by the final
+ // track
+ mutable int m_nb_track_hits; // number of track hits in the final
+ // track
+ mutable double m_chi2; // chi^2 of the final track
+ double m_m_x1, m_b_x1; // slope and intercept in the x1-x3 plane
+ double m_m_x2, m_b_x2; // slope and intercept in the x2-x3 plane
+ double m_m_x1_err, m_b_x1_err; // errors on m_x1 and b_x1
+ double m_m_x2_err, m_b_x2_err; // errors on m_x2 and b_x2
+ mutable MTStraightLine m_track; // final track
+
+// parameters for the adjustment of the track reconstruction //
+ double m_r_max; // maximum radius
+ double m_road_width; // road width for pattern recognition
+ double m_time_out; // time-out for track finding
+
+// chi^2 refitter //
+ DCSLFitter m_nfitter; // NIKHEF straight line reconstruction
+
+// internal storage vectors //
+ std::vector<MTStraightLine> m_tangent;
+ // vector of tangents (for 1 track)
+ std::vector<int> m_nb_cand_hits; // number of hits on track candidate
+ std::vector<MTStraightLine> m_candidate; // track candidates
+
+// initialization methods //
+ void m_init(void);
+ // default initialization: road width = 0.5 CLHEP::mm
+ void m_init(const double & r_road_width);
+ // initialization with user-defined road width
+
+// auxiliary methods //
+ MTStraightLine tangent(const Amg::Vector3D & r_w1,
+ const double & r_r1, const double & r_sigma12,
+ const Amg::Vector3D & r_w2, const double & r_r2,
+ const double & r_sigma22, const int & r_case) const;
+ // method for the calculation of tangents with errors;
+ // r_w1: wire position for the first hit,
+ // r_r1: drift radius of the first hit,
+ // r_sigma12: sigma(r_r1)^2,
+ // r_w2: wire position for the second hit,
+ // r_r2: drift radius of the second hit,
+ // r_sigma22: sigma(r_r2)^2,
+ // r_case = 1, 2, 3, 4: select one of the four cases of a tangent
+
+ MTStraightLine m_track_candidate(
+ const IndexSet & r_index_set,
+ const int & r_k_cand,
+ const int & r_l_cand,
+ const int & r_cand_case,
+ std::vector<Amg::Vector3D> r_w,
+ std::vector<double> r_r,
+ std::vector<double> r_sigma2,
+ double & r_chi2) const;
+ // method for the calculatation of a straight line from tangents;
+ // r_nb_selected_hits: number of selected hits,
+ // r_k_cand: index of the 1st wire of the candidate tangent
+ // (starting from 0),
+ // r_l_cand: index of the 2nd wire of the candidate tangent
+ // (starting from 0),
+ // r_cand_case: configuration of the tangent to be used (1, 2, 3, or 4)
+ // r_w: wire positions,
+ // r_r: drift radii,
+ // r_sigma2: sigma(r)^2,
+ // r_chi2: chi^2 for the reconstructed track
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/StraightPatRec.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/StraightPatRec.h
new file mode 100644
index 0000000000000000000000000000000000000000..9309b2460d4a1fc239628f55c1378ec7c8d55dec
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/MdtCalibFitters/StraightPatRec.h
@@ -0,0 +1,207 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 27.10.2007, AUTHOR: OLIVER KORTNER
+// Modified: 26.11.2007 by O. Kortner, fix for segment refinement.
+// 13.12.2007 by O. Kortner, time-out added.
+// 07.08.2008 by O. Kortner, bug fix when hits are disabled.
+// 19.11.2008 by I. Potrap: 1) several bugs have been fixed;
+// 2) external fitter(DCLS) replaced by
+// internal code of linear-least-chi2-fit
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+#ifndef StraightPatRecH
+#define StraightPatRecH
+
+//::::::::::::::::::::::::::
+//:: CLASS StraightPatRec ::
+//::::::::::::::::::::::::::
+
+/// \class StraightPatRec
+/// This class searches for the best hit pattern to be fitted with the
+/// DCSLFitter.
+/// The user can set a time-out for the track finding. It is set to 10 seconds
+/// by default.
+///
+/// \author Oliver.Kortner@CERN.CH
+///
+/// \date 27.10.2007, 13.12.2007
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+// CLHEP //
+//#include "CLHEP/config/CLHEP.h"
+//#include "CLHEP/Units/SystemOfUnits.h"
+//#include "CLHEP/Units/PhysicalConstants.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+
+// STL //
+#include <vector>
+
+// MuonCalib //
+#include "MdtCalibInterfaces/IMdtPatRecFitter.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "MdtCalibFitters/MTStraightLine.h"
+#include "MdtCalibFitters/DCSLFitter.h"
+//#include "MdtCalibFitters/FourCasesFitter.h"
+
+
+namespace MuonCalib {
+
+class StraightPatRec : public IMdtPatRecFitter {
+
+public:
+// Constructors //
+ StraightPatRec(void) {
+ m_init();
+ }
+ ///< Default constructor: road width for pattern recognition = 0.5 mm.
+
+ StraightPatRec(const double & r_road_width) {
+ m_init(r_road_width);
+ }
+ ///< Constructor: user-defined road width for pattern recognition.
+
+// Methods //
+// get-methods //
+ double roadWidth(void) const;
+ ///< get the road width used in the
+ ///< pattern recognition
+ unsigned int numberOfTrackHits(void) const;
+ ///< get the number of track hits in the
+ ///< final track,
+ ///< the method returns 0, if no track
+ ///< has been reconstructed
+ const std::vector<const MdtCalibHitBase*> & trackHits(void) const;
+ ///< get a vector with a pointer to the
+ ///< hits used in the track
+ ///< reconstruction
+ double chi2(void) const;
+ ///< get the chi^2 of the final track,
+ ///< the method returns -1, if no track
+ ///< has been reconstructed
+ double chi2PerDegreesOfFreedom(void) const;
+ ///< get the chi^2 per degrees of
+ ///< freedom for the final track,
+ ///< the method returns -1, if no track
+ ///< has been reconstructed
+ MTStraightLine track(void) const;
+ ///< get the final track in the local
+ ///< co-ordinate frame, i.e.
+ ///< z
+ ///< ^
+ ///< o o o o o o |
+ ///< ... o o o o o ... o--> y
+ ///< o o o o o o x
+
+// set-method //
+ void setRoadWidth(const double & r_road_width);
+ ///< set the road width for the pattern
+ ///< recognition = r_road_width
+ void setTimeOut(const double & time_out);
+ ///< set the time-out for the track
+ ///< finding to time_out (in seconds)
+
+ void setFixSelection(bool fix_sel);
+
+// methods required by the base class "IMdtSegmentFitter" //
+ bool fit(MuonCalibSegment & r_segment) const;
+ ///< reconstruction of the track using
+ ///< all hits in the segment
+ ///< "r_segment", returns true in case
+ ///< of success;
+ ///< the algorithm overwrites
+ ///< the track radii, the track
+ ///< position, track direction, and
+ ///< chi^2 per degrees of freedom;
+ ///< warning: the errors of the track
+ ///< radii are only approximate
+ bool fit(MuonCalibSegment & r_segment, HitSelection r_selection) const;
+ ///< reconstruction of the track using
+ ///< only those hits in r_segment for
+ ///< which the r_selection[.] is 0,
+ ///< return true in case of success;
+ ///< the algorithm overwrites
+ ///< the track position, direction,
+ ///< and the chi^2 in r_segment; it
+ ///< updates the distances of all hits
+ ///< from the track, i.e. also of those
+ ///< hits which were rejected from the
+ ///< track reconstruction;
+ ///< warning : the errors of the track
+ ///< CLHEP::radii are only approximate
+
+bool fitCallByReference(MuonCalibSegment & r_segment, HitSelection &r_selection) const;
+
+ void printLevel(int level);
+ ///< set the print level, this call
+ ///< has no effect
+
+private:
+// internal co-ordinate definition //
+// x3
+// ^
+// o o o o o o |
+// ... o o o o o ... o--> x2
+// o o o o o o x1
+//
+
+// final track parameters //
+ mutable std::vector<const MdtCalibHitBase*> m_track_hits;
+ // pointer to hits used by the final
+ // track
+ mutable int m_nb_track_hits; // number of track hits in the final
+ // track
+ mutable double m_chi2; // chi^2 of the final track
+ double m_m_x1, m_b_x1; // slope and intercept in the x1-x3 plane
+ double m_m_x2, m_b_x2; // slope and intercept in the x2-x3 plane
+ double m_m_x1_err, m_b_x1_err; // errors on m_x1 and b_x1
+ double m_m_x2_err, m_b_x2_err; // errors on m_x2 and b_x2
+ mutable MTStraightLine m_track; // final track
+
+ bool m_fix_selection;
+
+// parameters for the adjustment of the track reconstruction //
+ double m_r_max; // maximum radius
+ double m_road_width; // road width for pattern recognition
+ double m_time_out; // time-out for track finding
+
+// chi^2 refitter //
+// DCSLFitter m_nfitter; // NIKHEF straight line reconstruction
+// FourCasesFitter m_nfitter;
+
+// initialization methods //
+ void m_init(void);
+ // default initialization: road width = 0.5 CLHEP::mm
+ void m_init(const double & r_road_width);
+ // initialization with user-defined road width
+
+// auxiliary methods //
+ MTStraightLine tangent(const Amg::Vector3D & r_w1,
+ const double & r_r1, const double & r_sigma12,
+ const Amg::Vector3D & r_w2, const double & r_r2,
+ const double & r_sigma22, const int & r_case) const;
+ // method for the calculation of tangents with errors;
+ // r_w1: wire position for the first hit,
+ // r_r1: drift radius of the first hit,
+ // r_sigma12: sigma(r_r1)^2,
+ // r_w2: wire position for the second hit,
+ // r_r2: drift radius of the second hit,
+ // r_sigma22: sigma(r_r2)^2,
+ // r_case = 1, 2, 3, 4: select one of the four cases of a tangent
+ MTStraightLine fitCandidate(MuonCalibSegment & r_segment,
+ const std::vector<unsigned int> & r_selection,
+ const MTStraightLine & cand_line) const;
+// const std::vector<unsigned int> & index) const;
+ // refit the candidate "cand_line" using the hits specified by
+ // the r_selection and index vectors
+
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/cmt/requirements b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/cmt/requirements
new file mode 100644
index 0000000000000000000000000000000000000000..5d7be1e3ea044bdb63e6922fed6a291d8b648ffb
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/cmt/requirements
@@ -0,0 +1,27 @@
+package MdtCalibFitters
+
+author Niels van Eldik
+
+use AtlasPolicy AtlasPolicy-*
+
+use GeoPrimitives GeoPrimitives-* DetectorDescription
+use EventPrimitives EventPrimitives-* Event
+
+use MuonCalibEventBase MuonCalibEventBase-* MuonSpectrometer/MuonCalib
+use MdtCalibInterfaces MdtCalibInterfaces-* MuonSpectrometer/MuonCalib/MdtCalib
+use MuonCalibMath MuonCalibMath-* MuonSpectrometer/MuonCalib/MuonCalibUtils
+
+include_dirs "$(MdtCalibFitters_root)"
+
+macro MdtCalibFitters_linkopts "-L$(MdtCalibFitters_root)/lib -lMdtCalibFitters"
+
+library MdtCalibFitters ../src/*.cxx
+
+apply_pattern installed_library
+
+private
+
+use AtlasCLHEP AtlasCLHEP-* External
+
+#macro cppdebugflags ''
+#macro_remove componentshr_linkopts "-Wl,-s"
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/Doxyfile b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/Doxyfile
new file mode 100644
index 0000000000000000000000000000000000000000..6d20c82d4f36ae64b076aa75a3e43cf1c5b6a284
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/Doxyfile
@@ -0,0 +1,160 @@
+# Doxyfile 1.2.10
+
+#---------------------------------------------------------------------------
+# General configuration options
+#---------------------------------------------------------------------------
+PROJECT_NAME = MdtCalibFitters
+PROJECT_NUMBER = MdtCalibFitters-00-00-03
+OUTPUT_DIRECTORY = /afs/cern.ch/user/d/domizia/calibration4/InstallArea/doc/MdtCalibFitters
+OUTPUT_LANGUAGE = English
+EXTRACT_ALL = YES
+EXTRACT_PRIVATE = YES
+EXTRACT_STATIC = YES
+HIDE_UNDOC_MEMBERS = NO
+HIDE_UNDOC_CLASSES = NO
+BRIEF_MEMBER_DESC = YES
+REPEAT_BRIEF = YES
+ALWAYS_DETAILED_SEC = NO
+FULL_PATH_NAMES = NO
+STRIP_FROM_PATH =
+INTERNAL_DOCS = NO
+CLASS_DIAGRAMS = YES
+SOURCE_BROWSER = YES
+INLINE_SOURCES = YES
+STRIP_CODE_COMMENTS = YES
+CASE_SENSE_NAMES = YES
+SHORT_NAMES = NO
+HIDE_SCOPE_NAMES = NO
+VERBATIM_HEADERS = YES
+SHOW_INCLUDE_FILES = YES
+JAVADOC_AUTOBRIEF = YES
+INHERIT_DOCS = YES
+INLINE_INHERITED_MEMB = YES
+INLINE_INFO = YES
+SORT_MEMBER_DOCS = NO
+DISTRIBUTE_GROUP_DOC = NO
+TAB_SIZE = 8
+GENERATE_TODOLIST = YES
+GENERATE_TESTLIST = YES
+GENERATE_BUGLIST = YES
+ALIASES =
+ENABLED_SECTIONS =
+MAX_INITIALIZER_LINES = 30
+OPTIMIZE_OUTPUT_FOR_C = NO
+SHOW_USED_FILES = YES
+#---------------------------------------------------------------------------
+# configuration options related to warning and progress messages
+#---------------------------------------------------------------------------
+QUIET = NO
+WARNINGS = YES
+WARN_IF_UNDOCUMENTED = YES
+WARN_FORMAT =
+WARN_LOGFILE =
+#---------------------------------------------------------------------------
+# configuration options related to the input files
+#---------------------------------------------------------------------------
+INPUT = ../src ../MdtCalibFitters ../doc
+FILE_PATTERNS = *.cxx *.h *.py *.xml *.mk
+RECURSIVE = YES
+EXCLUDE =
+EXCLUDE_PATTERNS =
+EXAMPLE_PATH = ../doc ../cmt ../share
+EXAMPLE_PATTERNS = *.cxx *.html requirements *.py
+IMAGE_PATH =
+INPUT_FILTER =
+FILTER_SOURCE_FILES = NO
+#---------------------------------------------------------------------------
+# configuration options related to the alphabetical class index
+#---------------------------------------------------------------------------
+ALPHABETICAL_INDEX = NO
+COLS_IN_ALPHA_INDEX = 5
+IGNORE_PREFIX =
+#---------------------------------------------------------------------------
+# configuration options related to the HTML output
+#---------------------------------------------------------------------------
+GENERATE_HTML = YES
+HTML_OUTPUT = html
+HTML_HEADER =
+HTML_FOOTER =
+HTML_STYLESHEET =
+HTML_ALIGN_MEMBERS = YES
+GENERATE_HTMLHELP = NO
+GENERATE_CHI = NO
+BINARY_TOC = NO
+TOC_EXPAND = NO
+DISABLE_INDEX = NO
+ENUM_VALUES_PER_LINE = 4
+GENERATE_TREEVIEW = NO
+TREEVIEW_WIDTH = 250
+#---------------------------------------------------------------------------
+# configuration options related to the LaTeX output
+#---------------------------------------------------------------------------
+GENERATE_LATEX = NO
+LATEX_OUTPUT =
+COMPACT_LATEX = NO
+PAPER_TYPE = a4wide
+EXTRA_PACKAGES =
+LATEX_HEADER =
+PDF_HYPERLINKS = NO
+USE_PDFLATEX = NO
+LATEX_BATCHMODE = NO
+#---------------------------------------------------------------------------
+# configuration options related to the RTF output
+#---------------------------------------------------------------------------
+GENERATE_RTF = NO
+RTF_OUTPUT =
+COMPACT_RTF = NO
+RTF_HYPERLINKS = NO
+RTF_STYLESHEET_FILE =
+RTF_EXTENSIONS_FILE =
+#---------------------------------------------------------------------------
+# configuration options related to the man page output
+#---------------------------------------------------------------------------
+GENERATE_MAN = NO
+MAN_OUTPUT =
+MAN_EXTENSION =
+MAN_LINKS = NO
+#---------------------------------------------------------------------------
+# configuration options related to the XML output
+#---------------------------------------------------------------------------
+GENERATE_XML = NO
+#---------------------------------------------------------------------------
+# Configuration options related to the preprocessor
+#---------------------------------------------------------------------------
+ENABLE_PREPROCESSING = YES
+MACRO_EXPANSION = NO
+EXPAND_ONLY_PREDEF = NO
+SEARCH_INCLUDES = YES
+INCLUDE_PATH =
+INCLUDE_FILE_PATTERNS =
+PREDEFINED =
+EXPAND_AS_DEFINED =
+#---------------------------------------------------------------------------
+# Configuration::addtions related to external references
+#---------------------------------------------------------------------------
+TAGFILES = /afs/cern.ch/atlas/software/dist/10.2.0/InstallArea/doc/AtlasPolicy.tag /afs/cern.ch/user/d/domizia/calibration4//InstallArea/doc/MuonCalibEvent.tag /afs/cern.ch/user/d/domizia/calibration4//InstallArea/doc/MuonCalibEventBase.tag
+GENERATE_TAGFILE = /afs/cern.ch/user/d/domizia/calibration4/InstallArea/doc/MdtCalibFitters.tag
+ALLEXTERNALS = NO
+PERL_PATH =
+#---------------------------------------------------------------------------
+# Configuration options related to the dot tool
+#---------------------------------------------------------------------------
+HAVE_DOT = YES
+CLASS_GRAPH = YES
+COLLABORATION_GRAPH = YES
+TEMPLATE_RELATIONS = YES
+INCLUDE_GRAPH = YES
+INCLUDED_BY_GRAPH = YES
+GRAPHICAL_HIERARCHY = YES
+DOT_PATH =
+DOTFILE_DIRS =
+MAX_DOT_GRAPH_WIDTH = 1024
+MAX_DOT_GRAPH_HEIGHT = 1024
+GENERATE_LEGEND = YES
+DOT_CLEANUP = YES
+DOT_IMAGE_FORMAT = gif
+UML_LOOK = YES
+#---------------------------------------------------------------------------
+# Configuration::addtions related to the search engine
+#---------------------------------------------------------------------------
+SEARCHENGINE = NO
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/mainpage.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/mainpage.h
new file mode 100644
index 0000000000000000000000000000000000000000..874015e5d5b5c1eaefccdded6bd45a2e855227e7
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/mainpage.h
@@ -0,0 +1,35 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+@mainpage MdtCalibFitters Package
+@author Niels.Van.Eldik@cern.ch, Fabrizio.Petrucci@cern.ch
+
+@section MdtCalibFittersIntro Introduction
+Package containing local pattern recognition algorithms and fitters that
+work on segment level.
+
+@section MdtCalibFittersOverview Class Overview
+
+- MuonCalib::DCSLFitter: straight line fitter for drift circles
+- MuonCalib::LocalSegmentResolver: calculates the local position and direction of the segment using the hits on track
+- MuonCalib::LocalToPrecision: provides transformation between local and precision reference frames
+- MuonCalib::MuCCaFitter: interface to the straight line fitter for drift circles used by Calib
+- MuonCalib::MuCCaFitterImplementation: the straight line fitter for drift circles used by Calib
+
+@ref used_MdtCalibFitters
+
+@ref requirements_MdtCalibFitters
+*/
+
+/**
+@page used_MdtCalibFitters Used Packages
+@htmlinclude used_packages.html
+*/
+
+/**
+@page requirements_MdtCalibFitters Requirements
+@include requirements
+*/
+
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/used_packages.html b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/used_packages.html
new file mode 100644
index 0000000000000000000000000000000000000000..9c7000d00b1dbf1f73bc890c2aedf5274e834762
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/doc/used_packages.html
@@ -0,0 +1,7 @@
+<center><table border='1'><tr><td><center><i>Used packages</i></center></td></tr>
+<tr><td><a href='../../AtlasCLHEP/html'>External/AtlasCLHEP</a></td></tr>
+<tr><td><a href='../../AtlasPolicy/html'>AtlasPolicy</a></td></tr>
+<tr><td><a href='../../MdtCalibInterfaces/html'>MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibInterfaces</a></td></tr>
+<tr><td><a href='../../MuonCalibEvent/html'>MuonSpectrometer/MuonCalib/MuonCalibEvent</a></td></tr>
+<tr><td><a href='../../MuonCalibEventBase/html'>MuonSpectrometer/MuonCalib/MuonCalibEventBase</a></td></tr>
+</table></center>
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedCandidateFinder.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedCandidateFinder.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d1102d8809f8535684bee60c5a1353fd8d3924e0
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedCandidateFinder.cxx
@@ -0,0 +1,193 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 13.07.2008, AUTHOR: OLIVER KORTNER
+// Modified: 04.08.2008 by O. Kortner, estimated direction of incidence can be
+// set.
+// 07.08.2008 by O. Kortner, bug fig in the pattern recognition.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS CurvedCandidateFinder ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibFitters/CurvedCandidateFinder.h"
+
+//::::::::::::::::::::::::
+//:: NAMESPACE SETTINGS ::
+//::::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: CONSTRUCTOR ::
+//:::::::::::::::::
+
+CurvedCandidateFinder::CurvedCandidateFinder(
+ const std::vector<const MdtCalibHitBase *> & hits) {
+
+ m_hits = hits;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD getCandidates ::
+//::::::::::::::::::::::::::
+
+const std::vector<CurvedLine> & CurvedCandidateFinder::getCandidates(
+ const double & road_width) {
+
+ Amg::Vector3D est_dir(0.0, 0.0, 1.0);
+ return getCandidates(road_width, est_dir);
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD getCandidates ::
+//::::::::::::::::::::::::::
+
+const std::vector<CurvedLine> & CurvedCandidateFinder::getCandidates(
+ const double & road_width,
+ const Amg::Vector3D & est_dir) {
+
+/////////////////////////////////////////
+// RETURN IF THERE ARE NOT ENOUGH HITS //
+/////////////////////////////////////////
+
+ if (m_hits.size()<3) {
+ cerr << endl
+ << "Class CurvedCandidateFinder, method getCandidates: "
+ << "WARNING!\n"
+ << "Not enough hits to determine a parabola!\n";
+ }
+
+///////////////
+// VARIABLES //
+///////////////
+
+ const MdtCalibHitBase *hit[3]; // three hits defining the candidate line
+ double min_z, max_z, dist(0.0); // auxialiary variables to define the points
+ vector<Amg::Vector3D> points(3); // points defining the curved candidate line
+ int sign[3]; // auxiliary sign array
+ Amg::Vector3D null(0.0, 0.0, 0.0); // auxiliary 0 vector
+ Amg::Vector3D xhat(1.0, 0.0, 0.0); // auxiliary unit vector
+ Amg::Vector3D shift_vec(0.0, est_dir.z(), -est_dir.y());
+ shift_vec = shift_vec.unit();
+
+/////////////////////////////////////////////
+// FIND THREE HITS WITH LARGEST SEPARATION //
+/////////////////////////////////////////////
+
+ hit[0] = m_hits[0];
+ hit[1] = 0;
+ hit[2] = m_hits[0];
+
+ min_z = hit[0]->localPosition().z();
+ max_z = hit[2]->localPosition().z();
+
+// get the points with the smallest and largest local z coordinate //
+ for (unsigned int k=1; k<m_hits.size(); k++) {
+ if (m_hits[k]->localPosition().z()<min_z) {
+ min_z = m_hits[k]->localPosition().z();
+ hit[0] = m_hits[k];
+ }
+ if (m_hits[k]->localPosition().z()>max_z) {
+ max_z = m_hits[k]->localPosition().z();
+ hit[2] = m_hits[k];
+ }
+ }
+
+// find a third hit with large separation from these two points //
+ for (unsigned int k=0; k<m_hits.size(); k++) {
+ if (m_hits[k]!=hit[0] && m_hits[k]!=hit[2]) {
+ if (hit[1]==0) {
+ hit[1] = m_hits[k];
+ dist = (hit[2]->localPosition().z()-hit[1]->localPosition().z())
+ -(hit[1]->localPosition().z()-hit[0]->localPosition().z());
+ } else {
+ if (dist<(hit[2]->localPosition().z()-
+ m_hits[k]->localPosition().z())
+ -(m_hits[k]->localPosition().z()-
+ hit[0]->localPosition().z())) {
+ dist = (hit[2]->localPosition().z()-
+ m_hits[k]->localPosition().z())
+ -(m_hits[k]->localPosition().z()-
+ hit[0]->localPosition().z());
+ hit[1] = m_hits[k];
+ }
+ }
+ }
+ }
+
+////////////////////////////////////////////////////////////////////
+// CALCULATE CANDIDATE LINES AND COUNT THE NUMBER OF HITS ON THEM //
+////////////////////////////////////////////////////////////////////
+
+// clear candidate vector //
+ m_candidates.clear();
+
+// search for the candidates //
+ for (sign[0]=-1; sign[0]<2; sign[0]=sign[0]+2) {
+ for (sign[1]=-1; sign[1]<2; sign[1]=sign[1]+2) {
+ for (sign[2]=-1; sign[2]<2; sign[2]=sign[2]+2) {
+
+// get a candidate //
+ unsigned int nb_hits(0);
+ for (unsigned int l=0; l<3; l++) {
+// points[l] = Amg::Vector3D(hit[l]->localPosition().x(),
+// hit[l]->localPosition().y()+
+// sign[l]*hit[l]->driftRadius(),
+// hit[l]->localPosition().z());
+ points[l] = hit[l]->localPosition()+
+ sign[l]*hit[l]->driftRadius()*shift_vec;
+ }
+ CurvedLine cand_line(points);
+
+// refine the candidate //
+ for (unsigned int l=0; l<3; l++) {
+ MTStraightLine tangent(
+ cand_line.getTangent(hit[l]->localPosition().z()));
+ Amg::Vector3D delta(tangent.positionVector()-hit[l]->localPosition());
+ Amg::Vector3D aux_dir(delta+(delta.dot(tangent.directionVector().unit())*tangent.directionVector().unit()));
+ aux_dir = hit[l]->driftRadius()*aux_dir.unit();
+ Amg::Vector3D mem(points[l]);
+ points[l][1] = (hit[l]->localPosition().y()+aux_dir.y());
+ points[l][2] = (hit[l]->localPosition().z()+aux_dir.z());
+ }
+
+// count the number of hits on the line within the given road width //
+ for (unsigned int k=0; k<m_hits.size(); k++) {
+ MTStraightLine w(Amg::Vector3D(0.0, m_hits[k]->localPosition().y(),
+ m_hits[k]->localPosition().z()), xhat, null, null);
+ double d(fabs((cand_line.getTangent(m_hits[k]->localPosition().z()
+ )).signDistFrom(w)));
+ if (fabs(m_hits[k]->driftRadius()-d)<road_width) {
+ nb_hits++;
+ }
+ }
+
+// add the candidate line to list of candidates if there are enough hits //
+ if (nb_hits==m_hits.size()) {
+ m_candidates.push_back(cand_line);
+ }
+
+ }
+ }
+ }
+
+ return m_candidates;
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedLine.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedLine.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..915f33c6f4cd538964dd7b169b47036e014c8e6a
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedLine.cxx
@@ -0,0 +1,223 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 05.04.2008, AUTHOR: OLIVER KORTNER
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS CurvedLine ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibFitters/CurvedLine.h"
+#include "MuonCalibMath/BaseFunctionFitter.h"
+#include "MuonCalibMath/LegendrePolynomial.h"
+
+//::::::::::::::::::::::::
+//:: NAMESPACE SETTINGS ::
+//::::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: DEFAULT CONSTRUCTOR ::
+//:::::::::::::::::::::::::
+
+CurvedLine::CurvedLine(void) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ vector<Amg::Vector3D> points(3);
+ vector<Amg::Vector3D> errors(3);
+
+////////////////////
+// FILL VARIABLES //
+////////////////////
+
+ points[0] = Amg::Vector3D(0.0, 0.0, 0.0); errors[0] = Amg::Vector3D(1.0, 1.0, 0.0);
+ points[1] = Amg::Vector3D(0.0, 0.0, 1.0); errors[1] = Amg::Vector3D(1.0, 1.0, 0.0);
+ points[2] = Amg::Vector3D(0.0, 0.0, 2.0); errors[2] = Amg::Vector3D(1.0, 1.0, 0.0);
+
+////////////////////
+// INITIALIZATION //
+////////////////////
+
+ init(points, errors);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: CONSTRUCTOR ::
+//:::::::::::::::::
+
+CurvedLine::CurvedLine(std::vector<Amg::Vector3D> & points) {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ vector<Amg::Vector3D> errors(points.size());
+
+////////////////////
+// FILL VARIABLES //
+////////////////////
+
+ for (unsigned int k=0; k<errors.size(); k++) {
+ errors[k] = Amg::Vector3D(1.0, 1.0, 0.0);
+ }
+
+////////////////////
+// INITIALIZATION //
+////////////////////
+
+ init(points, errors);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: CONSTRUCTOR ::
+//:::::::::::::::::
+
+CurvedLine::CurvedLine(std::vector<Amg::Vector3D> & points,
+ std::vector<Amg::Vector3D> & x_and_y_errors) {
+
+ init(points, x_and_y_errors);
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::
+//:: METHOD getPointOnLine ::
+//::::::::::::::::::::::::::::::
+
+Amg::Vector3D CurvedLine::getPointOnLine(const double & loc_z) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ double loc_x(0.0), loc_y(0.0);
+
+///////////////////////////////
+// CALCULATE THE COORDINATES //
+///////////////////////////////
+
+ for (int k=0; k<m_coeff_xz.rows(); k++) {
+ loc_x = loc_x+m_coeff_xz[k]*m_Legendre->value(k, loc_z);
+ }
+
+ for (int k=0; k< m_coeff_yz.rows(); k++) {
+ loc_y = loc_y+m_coeff_yz[k]*m_Legendre->value(k, loc_z);
+ }
+
+////////////////////////////////
+// RETURN THE REQUESTED POINT //
+////////////////////////////////
+
+ return Amg::Vector3D(loc_x, loc_y, loc_z);
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD getTangent ::
+//::::::::::::::::::::::::
+
+MTStraightLine CurvedLine::getTangent(const double & loc_z) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ Amg::Vector3D null_vec(0.0, 0.0, 0.0); // auxiliary 0 vector
+ Amg::Vector3D point_1(getPointOnLine(loc_z)); // first point on the curved
+ // line
+ Amg::Vector3D point_2(getPointOnLine(loc_z+1.0)); // second point on the
+ // curved line
+
+////////////////////////
+// RETURN THE TANGENT //
+////////////////////////
+
+ return MTStraightLine(point_1, point_2-point_1, null_vec, null_vec);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD init ::
+//:::::::::::::::::
+
+void CurvedLine::init(std::vector<Amg::Vector3D> & points,
+ std::vector<Amg::Vector3D> & x_and_y_errors) {
+
+////////////////////////////////
+// CHECK THE NUMBER OF POINTS //
+////////////////////////////////
+
+ if (points.size()<3) {
+ cerr << endl
+ << "Class CurvedLine, method init: ERROR!\n"
+ << "Not enough points given, must be at least 3 points!\n";
+ }
+
+///////////////
+// VARIABLES //
+///////////////
+
+ BaseFunctionFitter fitter;
+ LegendrePolynomial legendre; // Legendre polynomial needed by the base
+ // function fitter
+ vector<SamplePoint> sample_points(points.size()); // sample points needed
+ // by the base function
+ // fitter
+/////////////////////////////////////////////
+// FILL THE VARIABLES AND PERFORM THE FITS //
+/////////////////////////////////////////////
+
+// xz plane //
+ for (unsigned int k=0; k<points.size(); k++) {
+ sample_points[k].set_x1(points[k].z());
+ sample_points[k].set_x2(points[k].x());
+ sample_points[k].set_error(x_and_y_errors[k].x());
+ }
+ fitter.set_number_of_coefficients(2);
+ fitter.fit_parameters(sample_points, 1, sample_points.size(), &legendre);
+ m_coeff_xz = fitter.coefficients();
+
+// yz plane //
+ for (unsigned int k=0; k<points.size(); k++) {
+ sample_points[k].set_x1(points[k].z());
+ sample_points[k].set_x2(points[k].y());
+ sample_points[k].set_error(x_and_y_errors[k].y());
+ }
+ fitter.set_number_of_coefficients(3);
+ fitter.fit_parameters(sample_points, 1, sample_points.size(), &legendre);
+ m_coeff_yz = fitter.coefficients();
+
+//////////////////////////////////////////////
+// GET A POINTER TO THE LEGENDRE POLYNOMIAL //
+//////////////////////////////////////////////
+
+ m_Legendre = Legendre_polynomial::get_Legendre_polynomial();
+
+ return;
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedPatRec.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedPatRec.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0d8cf9eed85f80d74f567883008f86619ed51497
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/CurvedPatRec.cxx
@@ -0,0 +1,573 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 05.04.2008, AUTHOR: OLIVER KORTNER
+// Modified: 25.07.2008 by O. Kortner, improved pattern recognition by using
+// the class "CurvedCandidateFinder"
+// 04.08.2008 by O. Kortner, further improvements of the pattern
+// recognition for large incidence angles.
+// 07.08.2008 by O. Kortner, bug fig in the pattern recognition.
+// 18.08.2008 by O. Kortner, update of chi^2 and segment position
+// and direction added.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS CurvedPatRec ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibFitters/CurvedPatRec.h"
+#include "MdtCalibFitters/CurvedCandidateFinder.h"
+#include "MuonCalibMath/Combination.h"
+#include "time.h"
+#include "cmath"
+
+//::::::::::::::::::::::::
+//:: NAMESPACE SETTINGS ::
+//::::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: DEFAULT CONSTRUCTOR ::
+//:::::::::::::::::::::::::
+
+CurvedPatRec::CurvedPatRec(void) {
+
+ m_chi2 = -1.0;
+ m_road_width = 0.5;
+ m_track_hits.clear();
+ m_time_out = 10;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: CONSTRUCTOR ::
+//:::::::::::::::::
+
+CurvedPatRec::CurvedPatRec(const double & road_width) {
+
+ m_chi2 = -1.0;
+ m_road_width = road_width;
+ m_track_hits.clear();
+ m_time_out = 10;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD roadWidth ::
+//::::::::::::::::::::::
+
+double CurvedPatRec::roadWidth(void) const {
+
+ return m_road_width;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::
+//:: METHOD numberOfTrackHits ::
+//::::::::::::::::::::::::::::::
+
+unsigned int CurvedPatRec::numberOfTrackHits(void) const {
+
+ return m_track_hits.size();
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD trackHits ::
+//::::::::::::::::::::::
+
+const std::vector<const MdtCalibHitBase*> & CurvedPatRec::trackHits(
+ void) const {
+
+ return m_track_hits;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD chi2 ::
+//:::::::::::::::::
+
+double CurvedPatRec::chi2(void) const {
+
+ return m_chi2;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::
+//:: METHOD chi2PerDegreesOfFreedom ::
+//::::::::::::::::::::::::::::::::::::
+
+double CurvedPatRec::chi2PerDegreesOfFreedom(void) const {
+
+ return m_chi2/static_cast<double>(m_track_hits.size()-3);
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD curvedTrack ::
+//::::::::::::::::::::::::
+
+const CurvedLine & CurvedPatRec::curvedTrack(void) const {
+
+ return m_curved_track;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD setRoadWidth ::
+//:::::::::::::::::::::::::
+
+void CurvedPatRec::setRoadWidth(const double & r_road_width) {
+
+ m_road_width = r_road_width;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD setTimeOut ::
+//:::::::::::::::::::::::
+
+void CurvedPatRec::setTimeOut(const double & time_out) {
+
+ m_time_out = time_out;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD fit(.) ::
+//:::::::::::::::::::
+
+bool CurvedPatRec::fit(MuonCalibSegment & r_segment) const {
+
+// select all hits //
+ HitSelection selection(r_segment.mdtHitsOnTrack(), 0);
+// call the other fit function //
+ return fit(r_segment, selection);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD fit(.,.) ::
+//:::::::::::::::::::::
+
+bool CurvedPatRec::fit(MuonCalibSegment & r_segment,
+ HitSelection r_selection) const {
+///////////////
+// VARIABLES //
+///////////////
+
+ time_t start, end; // start and end times (needed for time-out)
+ double diff; // difference of start and end time (needed for time-out)
+ Combination combination;
+ vector<unsigned int> hit_index; // hit indices for a given combination
+ unsigned int try_nb_hits; // try to find a segment with try_nb_hits hits
+ bool segment_found(false); // flag indicating the a segment has been found
+ vector<const MdtCalibHitBase *> cand_track_hits; // vector of the track hits
+ // found so far
+ CurvedLine aux_line; // memory for reconstructed curved lines
+ Amg::Vector3D null(0.0, 0.0, 0.0);
+ Amg::Vector3D xhat(1.0, 0.0, 0.0);
+ vector<Amg::Vector3D> points; // hit points for the track fit
+ vector<const MdtCalibHitBase*> loc_track_hits; // track hit store
+
+////////////
+// RESETS //
+////////////
+
+ m_chi2 = -1.0;
+ m_track_hits.clear();
+ time(&start);
+
+////////////////////////////////////////
+// CHECK SIZE OF THE SELECTION VECTOR //
+////////////////////////////////////////
+
+ if (r_selection.size()!=r_segment.mdtHitsOnTrack()) {
+ cerr << endl
+ << "Class CurvedPatRec, method fit: ERROR!\n"
+ << "Size of selection vector does not match the number of hits on"
+ << "track!\n";
+ exit(1);
+ }
+
+//////////////////////
+// PREPARATORY WORK //
+//////////////////////
+
+// perform a straight track fit to get an estimate of the incidence angle //
+ Amg::Vector3D est_dir(0.0, 0.0, 1.0);
+ m_sfitter.setRoadWidth(2.0*m_road_width);
+ m_sfitter.setTimeOut(0.5*m_time_out);
+ if (m_sfitter.fit(r_segment, r_selection)) {
+ est_dir = m_sfitter.track().directionVector();
+ }
+
+// store track hits //
+ for (unsigned int k=0; k<r_segment.mdtHitsOnTrack(); k++) {
+ if (r_selection[k]==0 && r_segment.mdtHOT()[k]->sigmaDriftRadius()<100) {
+ m_track_hits.push_back(r_segment.mdtHOT()[k]);
+ loc_track_hits.push_back(r_segment.mdtHOT()[k]);
+ }
+ }
+
+// return, if there are too few hits //
+ if (m_track_hits.size()<4) {
+ return false;
+ }
+/////////////////////////
+// PATTERN RECOGNITION //
+/////////////////////////
+
+// try to find a segment with as many hits on it as possible //
+ try_nb_hits = loc_track_hits.size();
+ while (!segment_found && try_nb_hits>3) {
+
+ // reset //
+ m_chi2 = -1.0;
+
+ // loop over the combinations //
+ combination.setNewParameters(loc_track_hits.size(), try_nb_hits);
+ for (unsigned int cb=0; cb<combination.numberOfCombinations(); cb++) {
+
+ // time-out //
+ time (&end);
+ diff = difftime (end,start);
+ if (diff>m_time_out) {
+ cerr << endl
+ << "Class CurvedPatRec: "
+ << "time-out for track finding after "
+ << m_time_out
+ << " seconds!\n";
+ return false;
+ }
+
+ // analyse the hit combination //
+ if (cb==0) {
+ combination.currentCombination(hit_index);
+ } else {
+ combination.nextCombination(hit_index);
+ }
+ vector<const MdtCalibHitBase *> track_hits;
+ for (unsigned int k=0; k<try_nb_hits; k++) {
+ track_hits.push_back(loc_track_hits[hit_index[k]-1]);
+ }
+
+ // find candidates //
+ CurvedCandidateFinder finder(track_hits);
+ const vector<CurvedLine> &candidates(finder.getCandidates(
+ m_road_width, est_dir));
+ if (candidates.size()==0) {
+ continue;
+ }
+
+ segment_found = true;
+
+ // store the track hits //
+ // m_track_hits = track_hits;
+
+ for (unsigned int cand=0; cand<candidates.size(); cand++) {
+ vector<Amg::Vector3D> errors(track_hits.size());
+ for (unsigned int k=0; k<errors.size(); k++) {
+ if (track_hits[k]->sigmaDriftRadius()>0.0) {
+ errors[k] = Amg::Vector3D(1.0,
+ track_hits[k]->sigmaDriftRadius(), 0.0);
+ } else {
+ errors[k] = Amg::Vector3D(1.0, 1.0, 0.0);
+ }
+ }
+
+ // get hit points //
+ points = getHitPoints(track_hits, candidates[cand]);
+
+ // fit a curved line through the points //
+ aux_line = CurvedLine(points, errors);
+
+ // calculate chi^2 //
+ double tmp_chi2(0.0);
+ for (unsigned int k=0; k<track_hits.size(); k++) {
+ MTStraightLine tang(m_curved_track.getTangent(
+ (track_hits[k]->localPosition()).z()));
+ MTStraightLine wire(Amg::Vector3D(0.0,
+ track_hits[k]->localPosition().y(),
+ track_hits[k]->localPosition().z()),
+ xhat, null, null);
+ double d(fabs(tang.signDistFrom(wire)));
+ if (track_hits[k]->sigma2DriftRadius()!=0) {
+ tmp_chi2 = tmp_chi2+
+ std::pow(d-track_hits[k]->driftRadius(), 2)/
+ track_hits[k]->sigma2DriftRadius();
+ } else {
+ tmp_chi2 = tmp_chi2+
+ std::pow(d-track_hits[k]->driftRadius(), 2)/0.01;
+ }
+ }
+
+ // compare chi^2 with chi^2 values found so far //
+ if (m_chi2<0) {
+ m_chi2 = tmp_chi2;
+ m_curved_track = aux_line;
+ // store the track hits //
+ m_track_hits = track_hits;
+ } else {
+ if (tmp_chi2<m_chi2) {
+ m_chi2 = tmp_chi2;
+ m_curved_track = aux_line;
+ m_track_hits = track_hits;
+ }
+ }
+
+ }
+
+ }
+
+ try_nb_hits = try_nb_hits-1;
+
+ }
+
+ if (segment_found==false) {
+ return false;
+ }
+
+///////////////////////////////
+// SECOND REFINED CURVED FIT //
+///////////////////////////////
+
+// get hit points //
+ points = getHitPoints(m_track_hits, m_curved_track);
+ vector<Amg::Vector3D> errors(m_track_hits.size());
+ for (unsigned int k=0; k<errors.size(); k++) {
+ if (m_track_hits[k]->sigmaDriftRadius()>0.0) {
+ errors[k] = Amg::Vector3D(1.0,
+ m_track_hits[k]->sigmaDriftRadius(), 0.0);
+ } else {
+ errors[k] = Amg::Vector3D(1.0, 1.0, 0.0);
+ }
+ }
+
+// fit a curved line through the points //
+ m_curved_track = CurvedLine(points, errors);
+
+/////////////////////
+// CALCULATE CHI^2 //
+/////////////////////
+
+ m_chi2 = 0.0;
+ for (unsigned int k=0; k<m_track_hits.size(); k++) {
+ MTStraightLine tang(m_curved_track.getTangent(
+ (m_track_hits[k]->localPosition()).z()));
+ MTStraightLine wire(Amg::Vector3D(0.0,
+ m_track_hits[k]->localPosition().y(),
+ m_track_hits[k]->localPosition().z()),
+ xhat, null, null);
+ double d(fabs(tang.signDistFrom(wire)));
+ if (m_track_hits[k]->sigma2DriftRadius()!=0) {
+ m_chi2 = m_chi2+ std::pow(d-m_track_hits[k]->driftRadius(), 2)/
+ m_track_hits[k]->sigma2DriftRadius();
+ } else {
+ m_chi2 = m_chi2+ std::pow(d-m_track_hits[k]->driftRadius(), 2)/0.01;
+ }
+ }
+
+//////////////////////////
+// UPDATE HIT RESIDUALS //
+//////////////////////////
+
+ MuonCalibSegment::MdtHitIt it = r_segment.mdtHOTBegin();
+ while(it!=r_segment.mdtHOTEnd()){
+
+ MdtCalibHitBase& hit = const_cast< MdtCalibHitBase& >( **it );
+
+ Amg::Vector3D pos(0.0, (hit.localPosition()).y(),
+ (hit.localPosition()).z());
+ MTStraightLine aux_line(pos, xhat, null, null);
+
+ MTStraightLine tang(m_curved_track.getTangent(pos.z()));
+
+ double dist(tang.signDistFrom(aux_line)); // track distance
+ double dist_err(1.0); // unknown error of the track distance
+ hit.setDistanceToTrack(dist, dist_err);
+
+ ++it;
+
+ }
+
+ if (std::isnan(m_chi2)) {
+ m_chi2=1.0e6;
+ }
+
+///////////////////////////////////////////////
+// UPDATE SEGMENT POSITION, DIRECTION, CHI^2 //
+///////////////////////////////////////////////
+
+ MTStraightLine tangent(m_curved_track.getTangent(
+ (r_segment.position()).z()));
+ r_segment.set(m_chi2/static_cast<double>(m_track_hits.size()-3),
+ tangent.positionVector(),
+ tangent.directionVector());
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD printLevel ::
+//:::::::::::::::::::::::
+
+void CurvedPatRec::printLevel(int level) {
+
+ cerr << "Class CurvedPatRec: method printLevel: " << level
+ << " has no effect!\n";
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD getHitPoint ::
+//:::::::::::::::::::::::::
+
+Amg::Vector3D CurvedPatRec::getHitPoint(const MdtCalibHitBase * hit,
+ const MTStraightLine & straight_track) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+// double dy, dz; // wire coordinates in the yz (precision) plane
+// double my, by; // slope and intercept of the given line in the yz
+ // (precision) plane
+// double d0; // distance of the track from the wire
+
+////////////////////
+// FILL VARIABLES //
+////////////////////
+
+// dy = (hit->localPosition()).y();
+// dz = (hit->localPosition()).z();
+// my = straight_track.m_x2();
+// by = straight_track.b_x2();
+
+///////////////////////////////////////////
+// CALCULATE THE POINT OF CLOSE APPROACH //
+///////////////////////////////////////////
+
+// double z0((dz-(by-dy)*my)/(1+my*my));
+// double y0(my*z0+by);
+
+/////////////////////////
+// CALCULATE HIT POINT //
+/////////////////////////
+
+// d0 = sqrt((y0-dy)*(y0-dy)+(z0-dz)*(z0-dz));
+
+
+
+ Amg::Vector3D point = straight_track.positionVector() + (straight_track.directionVector().unit().dot(hit->localPosition() - straight_track.positionVector() ) ) * straight_track.directionVector().unit();
+ Amg::Vector3D point_2 = hit->localPosition() + hit->driftRadius() * (point - hit->localPosition()).unit() ;
+
+ return point_2;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD getHitPoints ::
+//:::::::::::::::::::::::::
+
+std::vector<Amg::Vector3D> CurvedPatRec::getHitPoints(
+ std::vector<const MdtCalibHitBase*> track_hits,
+ const MTStraightLine & straight_track) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ vector<Amg::Vector3D> hit_vec;
+
+/////////////////////
+// FILL HIT VECTOR //
+/////////////////////
+
+ for (unsigned int k=0; k<track_hits.size(); k++) {
+ hit_vec.push_back(getHitPoint(track_hits[k], straight_track));
+ }
+
+///////////////////////////
+// RETURN THE HIT VECTOR //
+///////////////////////////
+
+ return hit_vec;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD getHitPoints ::
+//:::::::::::::::::::::::::
+
+std::vector<Amg::Vector3D> CurvedPatRec::getHitPoints(
+ std::vector<const MdtCalibHitBase*> track_hits,
+ const CurvedLine & curved_track) const {
+///////////////
+// VARIABLES //
+///////////////
+
+ vector<Amg::Vector3D> hit_vec;
+
+/////////////////////
+// FILL HIT VECTOR //
+/////////////////////
+
+ for (unsigned int k=0; k<track_hits.size(); k++) {
+ hit_vec.push_back(getHitPoint(track_hits[k],
+ curved_track.getTangent(
+ (track_hits[k]->localPosition()).z())));
+ }
+
+///////////////////////////
+// RETURN THE HIT VECTOR //
+///////////////////////////
+
+ return hit_vec;
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..35e6c621045ded14b6f83103442d24ec8c1807a6
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/DCSLFitter.cxx
@@ -0,0 +1,277 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MdtCalibFitters/DCSLFitter.h"
+#include "cmath"
+
+namespace MuonCalib {
+
+ void DCSLFitter::printLevel(int level) {
+ if(level > 0 ) m_debug = true;
+ }
+
+ bool DCSLFitter::fit( MuonCalibSegment& seg ) const
+ {
+ // select all hits
+ HitSelection selection(seg.mdtHitsOnTrack(),0);
+
+ // call fit function
+ return fit(seg,selection);
+ }
+
+
+ bool DCSLFitter::fit( MuonCalibSegment& seg, HitSelection selection ) const
+ {
+ if(m_debug) std::cout << "New seg: " << std::endl; //<< seg;
+
+
+ int N = seg.mdtHitsOnTrack();
+
+ if(N<2) {
+ return false;
+ }
+
+ if((int)selection.size() != N){
+ selection.clear();
+ selection.assign(N,0);
+ }else{
+ int used(0);
+ for(int i=0;i<N;++i){
+ if(selection[i] == 0) ++used;
+ }
+ if(used < 2){
+ std::cout << "TO FEW HITS SELECTED" << std::endl;
+ return false;
+ }
+ }
+
+ Amg::Vector3D pos = seg.position();
+ Amg::Vector3D dir = seg.direction();
+
+ double S(0),Sz(0),Sy(0);
+ double Zc(0),Yc(0);
+
+ std::vector<double> y(N);
+ std::vector<double> z(N);
+ std::vector<double> r(N);
+ std::vector<double> w(N);
+ {
+ MuonCalibSegment::MdtHitIt hit = seg.mdtHOTBegin();
+
+ int ii(0);
+ while(hit != seg.mdtHOTEnd()){
+ const MdtCalibHitBase& h = **hit;
+
+ y[ii] = getY( h.localPosition() );
+ z[ii] = getZ( h.localPosition() );
+ r[ii] = std::abs( h.driftRadius() );
+ if(h.sigma2DriftRadius() > 0)
+ w[ii] = 1./( h.sigma2DriftRadius() );
+ else
+ w[ii] = 0.;
+ if(r[ii]<0){
+ r[ii] = 0.;
+ std::cout << "DCSLFitter ERROR: <Negative r> " << r[ii] << std::endl;
+ }
+ if(m_debug)
+ std::cout << "DC: (" << y[ii] << "," << z[ii] << ") R = " << r[ii]
+ << " W " << w[ii] << std::endl;
+
+ if(selection[ii]){
+ ++hit;++ii;
+ continue;
+ }
+ S+=w[ii];
+ Sz+= w[ii]*z[ii];
+ Sy+= w[ii]*y[ii];
+ ++hit;++ii;
+ }
+ }
+ Zc = Sz/S;
+ Yc = Sy/S;
+
+ if(m_debug)
+ std::cout << "Yc " << Yc << " Zc " << Zc << std::endl;
+ //
+ // shift hits
+ //
+ //std::cout << "Shifted d " << d << std::endl;
+ Sy = 0;
+ Sz = 0;
+ double Szz(0),Syy(0),Szy(0),Syyzz(0);
+ std::vector<double> rw(N);
+ std::vector<double> ryw(N);
+ std::vector<double> rzw(N);
+
+ for(int i=0;i<N;++i){
+
+ y[i] -= Yc;
+ z[i] -= Zc;
+
+ if(selection[i]) continue;
+
+ rw[i] = r[i]*w[i];
+ ryw[i] = rw[i]*y[i];
+ rzw[i] = rw[i]*z[i];
+
+ Szz += z[i]*z[i]*w[i];
+ Syy += y[i]*y[i]*w[i];
+ Szy += y[i]*z[i]*w[i];
+ Syyzz += (y[i]-z[i])*(y[i]+z[i])*w[i];
+ }
+
+ if(m_debug)
+ std::cout << " Szz " << Szz << " Syy " << Syy
+ << " Szy " << Szy << " Syyzz " << Syyzz << std::endl;
+
+ int count(0);
+ double R(0),Ry(0),Rz(0);
+ double Att(0);
+ double Add = S;
+ double Bt(0);
+ double Bd(0);
+ double Stt(0);
+ double Sd(0);
+
+ double cosin = getZ( dir )/dir.mag();
+ double sinus = getY( dir )/dir.mag();
+
+ if(m_debug)
+ std::cout << "cos " << cosin << " sin " << sinus << std::endl;
+
+ // make sure 0 <= theta < PI
+ if ( sinus < 0.0 ) {
+ if(m_debug) {
+ std::cout << "sinus < 0.0 " << std::endl;
+ }
+ sinus = -sinus;
+ cosin = -cosin;
+ } else if ( sinus == 0.0 && cosin < 0.0 ) {
+ if(m_debug) {
+ std::cout << "sinus == 0.0 && cosin < 0.0" << std::endl;
+ }
+ cosin = -cosin;
+ }
+ //
+ // calculate shift
+ //
+ double d = -( getZ( pos )-Zc)*sinus+( getY( pos )-Yc)*cosin;
+ double theta = acos(cosin);
+ if(m_debug) {
+ std::cout << "____________INITIAL VALUES________________" << count << std::endl;
+ std::cout << "Theta " << theta << " d " << d << std::endl;
+ }
+
+ while(count<100){
+ if(m_debug) std::cout << "____________NEW ITERATION________________" << count << std::endl;
+ R=0;Ry=0;Rz=0;
+ for(int i=0;i<N;++i){
+ if(selection[i]) continue;
+
+ double dist = y[i]*cosin-z[i]*sinus;
+ if(dist>d){
+ R -= rw[i];
+ Ry -= ryw[i];
+ Rz -= rzw[i];
+ if(m_debug) std::cout << " < - > " << dist - d << " r " << r[i];
+ }else{
+ R += rw[i];
+ Ry += ryw[i];
+ Rz += rzw[i];
+ if(m_debug) std::cout << " < + > " << dist - d << " r " << r[i];
+ }
+ }
+ if(m_debug) std::cout << std::endl;
+ Att = Syy + cosin*(2*sinus*Szy - cosin*Syyzz);
+ Bt = -Szy + cosin*(sinus*Syyzz + 2*cosin*Szy + Rz) + sinus*Ry;
+ Bd = -S*d + R;
+ if(Att==0){
+ std::cerr << "===> Error NewtonSLDCFitter ZERO Determinant" << std::endl;
+ return false;
+ }
+ theta += Bt/Att;
+ if ( theta < 0.0 ) theta += M_PI;
+ if ( theta >= M_PI ) theta -= M_PI;
+ cosin = cos(theta);
+ sinus = sqrt(1-cosin*cosin);
+ d = R/S;
+ if(m_debug) std::cout << "R " << R << " Ry " << Ry << " Rz " << Rz << std::endl;
+ if(m_debug) std::cout << "Att " << Att << " Add " << Add << " Bt " << Bt << " Bd " << Bd << std::endl;
+ if(m_debug) std::cout << "dTheta " << Bt/Att << " dD " << Bd/Add << std::endl;
+ if(fabs(Bt/Att) < 0.001 && fabs(Bd/Add) < 0.001){
+ Stt = sqrt(1/Att);
+ Sd = sqrt(1/Add);
+ if(m_debug) std::cout << "Fit converged after " << count << " iterations " << std::endl;
+
+ if(m_debug) std::cout << "Theta " << theta << "d " << d << std::endl;
+ if(m_debug) std::cout << "Errors: theta " << Stt << " d " << Sd << std::endl;
+
+// error_dy0 = Sd;
+// error_dtheta = Stt;
+ seg.setErrors(Sd,Stt);
+
+ break;
+ }
+ ++count;
+ }
+ if(m_debug) std::cout << "Calculating chi2" << std::endl;
+
+ std::vector< double > yl(N);
+ std::vector< double > dyl(N);
+ std::vector< double > res(N);
+ double chi2(0);
+ if(m_debug) std::cout << "contributions to chi2: " << std::endl;
+
+ // calculate predicted hit positions from track parameters
+ for(int i=0;i<N;++i){
+ yl[i] = cosin*y[i] - sinus*z[i] - d;
+ double dth = -(sinus*y[i] + cosin*z[i])*Stt;
+ dyl[i] = sqrt( dth*dth + Sd*Sd );
+ res[i] = fabs(yl[i]) - r[i];
+ if(m_debug) std::cout << " r_track " << yl[i] << " dr " << dyl[i]
+ << " r_rt " << r[i] << " res " << res[i]
+ << " pull " << res[i]*res[i]*w[i] << std::endl;
+ //skip hits that are not used in fit
+ if(selection[i]) continue;
+
+ chi2 += res[i]*res[i]*w[i];
+
+ }
+
+ // std::cout << " Chi2 = " << chi2 << " chi2 per dof " << chi2/(N-2) << std::endl;
+
+ if(m_debug) std::cout << "Fit complete: Chi2 tof " << chi2/(N-2)
+ << " " << !(chi2/(N-2) > 5) << std::endl;
+ if(chi2/(N-2) > 5) {
+ if(m_debug) {
+ std::cout << "_______NOT GOOD " << std::endl;
+ }
+ }
+
+ if(m_debug) std::cout << "Transforming back to real world" << std::endl;
+
+ Amg::Vector3D ndir = getVec( dir.x() , sinus, cosin );
+ Amg::Vector3D npos = getVec( pos.x() , Yc + cosin*d, Zc - sinus*d );
+
+ if(m_debug) std::cout << "New line: position " << npos << " direction " << ndir << std::endl;
+
+ seg.set( chi2/(N-2), npos, ndir );
+
+ MuonCalibSegment::MdtHitIt it = seg.mdtHOTBegin();
+
+ int i(0);
+ while(it!=seg.mdtHOTEnd()){
+ MdtCalibHitBase& hit = const_cast< MdtCalibHitBase& >( **it );
+
+ hit.setDistanceToTrack( yl[i], dyl[i] );
+
+ ++it;++i;
+ }
+
+ if(m_debug) std::cout << "fit done" << std::endl;
+ //tracking failed if position and directins are not real numbers
+ return (std::isfinite(ndir.y()) && std::isfinite(ndir.z()) && std::isfinite(npos.y()) && std::isfinite(npos.z()));
+ }
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/IndexSet.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/IndexSet.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..1f088ecfa42baafcd945f69385cf2dd8f3350d42
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/IndexSet.cxx
@@ -0,0 +1,211 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 09.05.2005, AUTHOR: OLIVER KORTNER
+// Modified: 16.07.2006 by O. Kortner, vector size check in initialization
+// added,
+// != operator corrected.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATIONS OF METHODS DEFINED IN THE CLASS IndexSet ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibFitters/IndexSet.h"
+#include <algorithm>
+#include <cstdlib>
+
+//:::::::::::::::::::::::
+//:: NAMESPACE SETTING ::
+//:::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD m_init ::
+//:::::::::::::::::::
+
+void IndexSet::m_init(void) {
+
+ m_nb_indices = 0;
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD m_init(.) ::
+//::::::::::::::::::::::
+
+void IndexSet::m_init(const unsigned int & r_nb_indices) {
+
+ m_nb_indices = r_nb_indices;
+ m_index = vector<int>(m_nb_indices);
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD m_init(.,.) ::
+//::::::::::::::::::::::::
+
+void IndexSet::m_init(const unsigned int & r_nb_indices,
+ const vector<int> r_index) {
+
+///////////////////////
+// CHECK VECTOR SIZE //
+///////////////////////
+
+ if (r_index.size()<r_nb_indices) {
+ cerr << endl
+ << "Class IndexSet, method m_init: ERROR!\n"
+ << "Index vector too short!\n";
+ exit(1);
+ }
+
+////////////////////
+// INITIALIZATION //
+////////////////////
+
+ m_nb_indices = r_nb_indices;
+ m_index = vector<int>(m_nb_indices);
+ for (unsigned int k=0; k<m_nb_indices; k++) {
+ if (k<r_index.size()) {
+ m_index[k] = r_index[k];
+ } else {
+ m_index[k] = 0;
+ }
+ }
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD size ::
+//:::::::::::::::::
+
+unsigned int IndexSet::size(void) const {
+
+ return m_nb_indices;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD resize ::
+//:::::::::::::::::::
+
+void IndexSet::resize(const unsigned int & r_size) {
+
+//:::::::::::::::
+//:: VARIABLES ::
+//:::::::::::::::
+
+ vector<int> aux_index = m_index; // vector for temporary storage of
+ // the indices
+
+//:::::::::::::::::::::::::::::::::::::
+//:: RESIZE AND COPY THE OLD INDICES ::
+//:::::::::::::::::::::::::::::::::::::
+
+ m_index = vector<int>(r_size);
+ for (unsigned int k=0; k<r_size; k++) {
+ if (k<m_nb_indices) {
+ m_index[k] = aux_index[k];
+ } else {
+ m_index[k] = 0;
+ }
+ }
+ m_nb_indices = r_size;
+
+ return;
+
+}
+
+//*****************************************************************************
+
+///:::::::::::::::::::::::
+//:: METHOD operator [] ::
+//::::::::::::::::::::::::
+
+int IndexSet::operator [] (const unsigned int & r_k) {
+
+ return m_index[r_k];
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD operator [] ::
+//::::::::::::::::::::::::
+
+int IndexSet::operator [] (const unsigned int & r_k) const {
+
+ return m_index[r_k];
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD sort ::
+//:::::::::::::::::
+
+void IndexSet::sort(void) {
+
+ std::sort(m_index.begin(), m_index.end());
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD operator == ::
+//::::::::::::::::::::::::
+
+bool IndexSet::operator == (const IndexSet & r_index_set) const {
+
+ if (r_index_set.size()!=m_nb_indices) {
+ return false;
+ }
+
+ for (unsigned int k=0; k<m_nb_indices; k++) {
+ if (m_index[k]!=r_index_set[k]) {
+ return false;
+ }
+ }
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::
+//:: METHOD != ::
+//:::::::::::::::
+
+bool IndexSet::operator != (const IndexSet & r_index_set) const {
+
+ return !(*this==r_index_set);
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/LocalSegmentResolver.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/LocalSegmentResolver.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..0c713c694b453f4e0c8c8f6af2514fcfefe5a75a
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/LocalSegmentResolver.cxx
@@ -0,0 +1,238 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MdtCalibFitters/LocalSegmentResolver.h"
+
+#include "MuonCalibEventBase/MdtCalibHitBase.h"
+
+#include "MdtCalibFitters/LocalToPrecision.h"
+
+#include <iostream>
+
+namespace MuonCalib {
+
+ LocalSegmentResolver::LocalSegmentResolver() : m_printLevel(0) {}
+
+ bool LocalSegmentResolver::resolve(MuonCalibSegment* seg) const{
+ if( !seg ){
+ if( m_printLevel >= 0 )
+ std::cout << "LocalSegmentResolver::resolve ERROR <got null pointer> " << std::endl;
+ return false;
+ }
+ if( seg->mdtHitsOnTrack() < 2 ) {
+ if( m_printLevel >= 1 )
+ std::cout << "LocalSegmentResolver::resolve ERROR <to few hits, cannot resolve direction> " << std::endl;
+ return false;
+ }
+
+ // get tangent lines
+ LineVec lines = getLines( *seg->mdtHOT().front(), *seg->mdtHOT().back() );
+
+ // find line which yields smallest residuals
+ int lineIndex = bestLine( seg->mdtHOT(), lines );
+
+ // check if everything went right
+ if( lineIndex < 0 ) return false;
+
+ // reset local position and direction of segment
+ seg->set( seg->chi2(), MuonCalib::LocalToPrecision::local(lines[lineIndex].first),
+ MuonCalib::LocalToPrecision::local(lines[lineIndex].second) );
+
+ return true;
+ }
+
+ LocalSegmentResolver::LineVec LocalSegmentResolver::getLines( const MdtCalibHitBase& firstHit,
+ const MdtCalibHitBase& lastHit ) const
+ {
+
+ Amg::Vector3D hpos1 = MuonCalib::LocalToPrecision::precision(firstHit.localPosition());
+ double x1 = hpos1.x();
+ double y1 = hpos1.y();
+ double r1 = std::abs(firstHit.driftRadius());
+ Amg::Vector3D hpos2 = MuonCalib::LocalToPrecision::precision(lastHit.localPosition());
+ double x2 = hpos2.x();
+ double y2 = hpos2.y();
+ double r2 = std::abs(lastHit.driftRadius());
+
+ double DeltaX = x2 - x1;
+ double DeltaY = y2 - y1;
+ double DistanceOfCenters = sqrt(DeltaX*DeltaX + DeltaY*DeltaY);
+ double Alpha0 = atan2(DeltaY,DeltaX);
+
+ LineVec list_of_lines;
+
+ if( m_printLevel >= 3 ){
+ std::cout << " calculating Lines (" << x1 << "," << y1 << ") " << r1
+ << " (" << x2 << "," << y2 << ") " << r2 << std::endl;
+ std::cout << " general dir " << (hpos2-hpos1).unit()
+ << " calculated : " << Amg::Vector3D( cos(Alpha0), sin(Alpha0), 0. ) << std::endl;
+ }
+
+ // Case of 0 drift distances, only 1 line
+ if ( r1 == 0. && r2 == 0.) {
+ Amg::Vector3D pos = hpos1;
+ Amg::Vector3D dir( cos(Alpha0) ,sin(Alpha0), 0 );
+ if( m_printLevel >= 5 ) std::cout << " line pos " << pos << " dir " << dir << std::endl;
+ list_of_lines.push_back( std::make_pair(pos, dir) );
+ return list_of_lines;
+ }
+
+
+ // Here are the first 2 "inner" lines ....
+ double RSum = r1 + r2;
+ double Alpha1 = asin(RSum/DistanceOfCenters);
+
+ double line_phi = Alpha0 + Alpha1;
+
+ Amg::Vector3D pos1( x1 + r1*sin(line_phi), y1 - r1*cos(line_phi), 0. ) ;
+ Amg::Vector3D dir1( cos(line_phi), sin(line_phi), 0. );
+
+ if( m_printLevel >= 5 )
+ std::cout << " line pos " << pos1 << " dir " << dir1 << std::endl;
+
+ list_of_lines.push_back( std::make_pair(pos1, dir1) );
+
+ line_phi = Alpha0 - Alpha1;
+
+ Amg::Vector3D pos2( x1 - r1*sin(line_phi), y1 + r1*cos(line_phi), 0. );
+ Amg::Vector3D dir2( cos(line_phi), sin(line_phi), 0. );
+
+ if( m_printLevel >= 5 )
+ std::cout << " line pos " << pos2 << " dir " << dir2 << std::endl;
+
+ list_of_lines.push_back( std::make_pair(pos2, dir2) );
+
+ // Case where one of the drifts is 0 ==> Only two lines
+ if (r1 == 0. || r2 == 0.) return list_of_lines;
+
+ // ... and here are the other 2 "outer" lines
+ double DeltaR = fabs(r2 - r1);
+ double Alpha2 = asin(DeltaR/DistanceOfCenters);
+
+ if ( r1 < r2 ) {
+ if( m_printLevel >= 5 )
+ std::cout << " r1 < r2 " << std::endl;
+
+
+ line_phi = Alpha0 + Alpha2;
+
+ Amg::Vector3D pos3( x1 - r1*sin(line_phi), y1 + r1*cos(line_phi), 0. );
+ Amg::Vector3D dir3( cos(line_phi),sin(line_phi), 0. );
+
+ if( m_printLevel >= 5 )
+ std::cout << " line pos " << pos3 << " dir " << dir3 << std::endl;
+
+ list_of_lines.push_back( std::make_pair(pos3, dir3) );
+
+ line_phi = Alpha0 - Alpha2;
+
+ Amg::Vector3D pos4( x1 + r1*sin(line_phi), y1 - r1*cos(line_phi), 0. );
+ Amg::Vector3D dir4( cos(line_phi),sin(line_phi), 0. );
+
+ if( m_printLevel >= 5 )
+ std::cout << " line pos " << pos4 << " dir " << dir4 << std::endl;
+
+ list_of_lines.push_back( std::make_pair(pos4, dir4) );
+
+ } else {
+ if( m_printLevel >= 5 )
+ std::cout << " r1 > r2 " << std::endl;
+
+ line_phi = Alpha0 + Alpha2;
+
+ Amg::Vector3D pos3( x1 + r1*sin(line_phi), y1 - r1*cos(line_phi), 0. );
+ Amg::Vector3D dir3( cos(line_phi),sin(line_phi), 0. );
+
+ if( m_printLevel >= 5 )
+ std::cout << " line pos " << pos3 << " dir " << dir3 << std::endl;
+
+ list_of_lines.push_back( std::make_pair(pos3, dir3) );
+
+ line_phi = Alpha0 - Alpha2;
+
+ Amg::Vector3D pos4( x1 - r1*sin(line_phi), y1 + r1*cos(line_phi), 0. );
+ Amg::Vector3D dir4( cos(line_phi),sin(line_phi), 0. );
+
+ if( m_printLevel >= 5 )
+ std::cout << " line pos " << pos4 << " dir " << dir4 << std::endl;
+
+ list_of_lines.push_back( std::make_pair(pos4, dir4) );
+ }
+
+ return list_of_lines;
+ }
+
+ int LocalSegmentResolver::bestLine( const LocalSegmentResolver::HitVec& hits,
+ const LocalSegmentResolver::LineVec& localTracks ) const
+ {
+ // loop over tangent lines
+ LineVec::const_iterator lit = localTracks.begin();
+ LineVec::const_iterator lit_end = localTracks.end();
+
+ // look for line with smallest residual sum
+ double ressummin = 1e20;
+ unsigned int resnum = 0;
+
+ for( ; lit!=lit_end;++lit ){
+
+ double ressum=0;
+
+ // calculate angle between line and z precision plane
+ double alpha = atan2(lit->second.y(),lit->second.x());
+
+ // get rotations in track frame
+ Amg::AngleAxis3D rotationAroundZ(-alpha,Amg::Vector3D::UnitZ());
+ //Amg::AngleAxis3D rotationAroundZ_inv(alpha,Amg::Vector3D::UnitZ());
+
+ // rotate position on track into track frame
+ Amg::Vector3D avePosTrk = rotationAroundZ*lit->first;
+
+ if( m_printLevel >= 5 ){
+ Amg::Vector3D lTrkDir = rotationAroundZ*lit->second;
+ std::cout << " angle " << alpha*57.32 << " trk dir in trk frame " << lTrkDir
+ << " pos " << avePosTrk << std::endl;
+ }
+ // loop over local hits
+ LocalSegmentResolver::HitVec::const_iterator it = hits.begin();
+ LocalSegmentResolver::HitVec::const_iterator it_end = hits.end();
+ for( ; it!=it_end; ++it ){
+
+ // rotate into track system
+ Amg::Vector3D spos = MuonCalib::LocalToPrecision::precision((*it)->localPosition());
+ Amg::Vector3D sposAve = rotationAroundZ*spos - avePosTrk;
+
+ // get drift radius (unsigned)
+ double r = std::abs( (*it)->driftRadius() );
+
+ // calculate residual and pull
+ double res = r - std::abs( sposAve.y() );
+
+ if( m_printLevel >= 5 )
+ std::cout << " r " << r << " r_trk " << std::abs( sposAve.y() ) << " residual " << res << std::endl;
+
+ ressum += res*res;
+ }
+
+ if( m_printLevel >= 3 )
+ std::cout << " line " << lit-localTracks.begin() << " residual sum " << ressum << std::endl;
+ if( ressum < ressummin ) {
+ ressummin = ressum;
+ resnum = lit-localTracks.begin();
+ }
+ }
+
+ if( m_printLevel >= 3 ){
+ std::cout << " Done selected line: ressum " << ressummin << " ## " << resnum << std::endl;
+ std::cout << " Position " << localTracks[resnum].first
+ << " direction " << localTracks[resnum].second << std::endl;
+ }
+ if( resnum >= localTracks.size() ) {
+ std::cout << " ERROR wrong line index " << std::endl;
+ return -1;
+ }
+
+ return resnum;
+ }
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/MTStraightLine.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/MTStraightLine.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a6a75483f943484aac6466af29cafa0a71ea478b
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/MTStraightLine.cxx
@@ -0,0 +1,333 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 01.03.2006, AUTHORS: OLIVER KORTNER, FELIX RAUSCHER
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibFitters/MTStraightLine.h"
+#include "cmath"
+#include <iostream>
+//:::::::::::::::::::::::
+//:: NAMESPACE SETTING ::
+//:::::::::::::::::::::::
+
+using namespace MuonCalib;
+
+//using namespace std;
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS MTStraightLine ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD m_init ::
+//:::::::::::::::::::
+
+void MTStraightLine::m_init(void) {
+
+ m_position = Amg::Vector3D(0.0, 0.0, 0.0);
+ m_direction = Amg::Vector3D(0.0, 0.0, 0.0);
+ m_position_error = Amg::Vector3D(0.0, 0.0, 0.0);
+ m_direction_error = Amg::Vector3D(0.0, 0.0, 0.0);
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD m_init(const Amg::Vector3D &, ...) ::
+//::::::::::::::::::::::::::::::::::::::::::::
+
+void MTStraightLine::m_init(const Amg::Vector3D & r_position,
+ const Amg::Vector3D & r_direction,
+ const Amg::Vector3D & r_position_error,
+ const Amg::Vector3D & r_direction_error) {
+
+ m_position = r_position;
+ m_direction = r_direction;
+ m_position_error = r_position_error;
+ m_direction_error = r_direction_error;
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD m_init ::
+//:::::::::::::::::::
+
+void MTStraightLine::m_init(const double & r_m_x1,
+ const double & r_b_x1,
+ const double & r_m_x2, const double & r_b_x2,
+ const double & r_m_x1_err, const double & r_b_x1_err,
+ const double & r_m_x2_err, const double & r_b_x2_err) {
+
+ m_position = Amg::Vector3D(r_b_x1, r_b_x2, 0.0);
+ m_direction = Amg::Vector3D(r_m_x1, r_m_x2, 1.0);
+ m_position_error = Amg::Vector3D(r_b_x1_err, r_b_x2_err, 0.0);
+ m_direction_error = Amg::Vector3D(r_m_x1_err, r_m_x2_err, 1.0);
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::
+//:: METHOD positionVector ::
+//:::::::::::::::::::::::::::
+
+Amg::Vector3D MTStraightLine::positionVector(void) const {
+
+ return m_position;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::
+//:: METHOD directionVector ::
+//::::::::::::::::::::::::::::
+
+Amg::Vector3D MTStraightLine::directionVector(void) const {
+
+ return m_direction;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::
+//:: METHOD positionError ::
+//::::::::::::::::::::::::::
+
+Amg::Vector3D MTStraightLine::positionError(void) const {
+
+ return m_position_error;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::
+//:: METHOD directionError ::
+//:::::::::::::::::::::::::::
+
+Amg::Vector3D MTStraightLine::directionError(void) const {
+
+ return m_direction_error;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD m_x1 ::
+//:::::::::::::::::
+
+double MTStraightLine::m_x1(void) const {
+
+ if (m_direction.z() == 0.0) {
+ return 0.0;
+ }
+ return m_direction.x()/m_direction.z();
+
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD m_x1_error ::
+//:::::::::::::::::::::::
+
+double MTStraightLine::m_x1_error(void) const {
+
+ if (m_direction.z()==0) {
+ return 0.0;
+ }
+
+ return sqrt(std::pow(m_direction_error.x()/m_direction_error.z(), 2)
+ +std::pow(m_direction_error.y()*m_x1()/m_direction_error.z(), 2));
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD b_x1 ::
+//:::::::::::::::::
+
+double MTStraightLine::b_x1(void) const {
+
+ if (m_direction.z() == 0.0) {
+ std::cerr << "Class MTStraightLine, method b_x1: WARNING!\n"
+ << "b_x1 not uniquely defined." << std::endl;
+ return m_position.x();
+ }
+ return (m_position.x()-m_position.z()*m_x1());
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD b_x1_error ::
+//:::::::::::::::::::::::
+
+double MTStraightLine::b_x1_error(void) const {
+
+ return sqrt(std::pow(m_position_error.x(), 2)+
+ std::pow(m_x1()*m_position_error.z(), 2)+
+ std::pow(m_position.z()*m_x1_error(),2));
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD m_x2 ::
+//:::::::::::::::::
+
+double MTStraightLine::m_x2(void) const {
+
+ if (m_direction.z() == 0.0) {
+ return 0.0;
+ }
+ return m_direction.y()/m_direction.z();
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD m_x2_error ::
+//:::::::::::::::::::::::
+
+double MTStraightLine::m_x2_error(void) const {
+
+ if (m_direction.z()==0) {
+ return 0.0;
+ }
+
+ return sqrt(std::pow(m_direction_error.y()/m_direction_error.z(), 2)
+ +std::pow(m_direction_error.y()*m_x2()/m_direction_error.z(), 2));
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD b_x2 ::
+//:::::::::::::::::
+
+double MTStraightLine::b_x2(void) const {
+
+ if (m_direction.z() == 0.0) {
+ std::cerr << "Class MTStraightLine, method b_x2: WARNING!\n"
+ << "b_x2 not uniquely defined." << std::endl;
+ return m_position.x();
+ }
+ return (m_position.y()-m_position.z()*m_x2());
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD b_x2_error ::
+//:::::::::::::::::::::::
+
+double MTStraightLine::b_x2_error(void) const {
+
+ return sqrt(std::pow(m_position_error.y(), 2)+
+ std::pow(m_x2()*m_position_error.z(), 2)+
+ std::pow(m_position.z()*m_x2_error(),2));
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::
+//:: METHOD pointOnLine ::
+//::::::::::::::::::::::::
+
+Amg::Vector3D MTStraightLine::pointOnLine(const double & lambda) const {
+
+ return m_position+lambda*m_direction;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD signDistFrom ::
+//:::::::::::::::::::::::::
+
+double MTStraightLine::signDistFrom(const MTStraightLine & h) const {
+
+//:::::::::::::::::::::::::
+//:: AUXILIARY VARIABLES ::
+//:::::::::::::::::::::::::
+
+ Amg::Vector3D a = m_position, u = m_direction;
+ // position and direction vectors of g
+ Amg::Vector3D b = h.positionVector(), v = h.directionVector();
+ // position and direction vectors of h
+ Amg::Vector3D n; // normal vector of plane spanned by g and h
+ Amg::Vector3D d; // distance vector
+
+//:::::::::::::
+//:: PROGRAM ::
+//:::::::::::::
+
+//::::::::::::::::::::::::
+//:: collinearity check ::
+//::::::::::::::::::::::::
+
+ n = u.cross(v);
+ d = a-b;
+ if (n.dot(n) == 0.0) {
+ // straight lines are parallel
+ // no sign given
+ return sqrt(d.dot(d)-(u.unit().dot(d))*(u.unit().dot(d)));
+ }
+
+ return (d.dot(n.unit()));
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD distFromLine ::
+//:::::::::::::::::::::::::
+
+double MTStraightLine::distFromLine(const Amg::Vector3D & point) const {
+
+//:::::::::::::::::::::::::
+//:: AUXILIARY VARIABLES ::
+//:::::::::::::::::::::::::
+
+ Amg::Vector3D u = m_direction;
+
+//::::::::::::::::::::::::
+//:: CALCULATE DISTANCE ::
+//::::::::::::::::::::::::
+
+ return sqrt((point-m_position).dot(point-m_position) -
+ ((point-m_position).dot(u.unit()))*((point-m_position).dot(u.unit())));
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/MuCCaFitter.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/MuCCaFitter.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..57a2b9b9721ebc995d150a1ed744a9a6b38a829c
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/MuCCaFitter.cxx
@@ -0,0 +1,653 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//#include <iostream.h>
+//#include <math.h>
+#include "MdtCalibFitters/MuCCaFitter.h"
+
+//class MuCCaFitterImplementation;
+
+namespace MuonCalib {
+
+void MuCCaFitter::printLevel(int level) {
+ if(level > 0 ) m_debug = true;
+}
+
+bool MuCCaFitter::fit( MuonCalibSegment& seg ) const
+{
+
+ // select all hits
+ HitSelection selection(seg.mdtHitsOnTrack(),0);
+
+ // call fit function
+ return fit(seg,selection);
+}
+
+
+bool MuCCaFitter::fit( MuonCalibSegment& seg, HitSelection selection ) const
+{
+ bool m_debug = false;
+ if(m_debug) std::cout << "New seg: " << std::endl; //<< seg;
+
+
+ int N = seg.mdtHitsOnTrack();
+
+ if(N<2) {
+ return false;
+ }
+
+ if((int)selection.size() != N){
+ selection.clear();
+ selection.assign(N,0);
+ }else{
+ int used(0);
+ for(int i=0;i<N;++i){
+ if(selection[i] == 0) ++used;
+ }
+ if(used < 2){
+ std::cout << "TO FEW HITS SELECTED" << std::endl;
+ return false;
+ }
+ }
+
+ Amg::Vector3D pos = seg.position();
+// Amg::Vector3D dir = seg.direction();
+
+ double S(0),Sy(0),Sz(0);
+ double Zc,Yc;
+ std::vector<double> y(N);
+ std::vector<double> x(N);
+ std::vector<double> z(N);
+ std::vector<double> r(N);
+ std::vector<double> sr(N);
+ std::vector<double> w(N);
+ std::vector<double> rw(N);
+ int ii;
+ int jj;
+ {
+ MuonCalibSegment::MdtHitIt hit = seg.mdtHOTBegin();
+
+ ii=0;
+ jj=0;
+ while(hit != seg.mdtHOTEnd()){
+ const MdtCalibHitBase& h = **hit;
+ y[ii] = getY( h.localPosition() );
+ z[ii] = getZ( h.localPosition() );
+ r[ii] = std::abs( h.driftRadius() );
+
+ if(h.sigma2DriftRadius() > 0){
+ w[ii] = 1./( h.sigma2DriftRadius() );
+ sr[ii] =sqrt( h.sigma2DriftRadius() );
+ } else {
+ w[ii] = 0.;
+ sr[ii]=0.;
+ }
+ if(m_debug)
+ std::cout << "MuCCaFitter: (" << x[ii] << "," << y[ii] << ") R = " << r[ii]
+ << " W " << w[ii] << std::endl;
+ rw[ii] = r[ii]*w[ii];
+ if(selection[jj]){
+ ++hit;++jj;
+ continue;
+ }
+ S+=w[ii];
+ Sz+= w[ii]*z[ii];
+ Sy+= w[ii]*y[ii];
+ ++hit;++ii;++jj;
+ }
+ }
+ Zc = Sz/S;
+ Yc = Sy/S;
+
+ MuCCaFitterImplementation *Fitter=new MuCCaFitterImplementation();
+ if(m_debug){
+ for (int i = 0; i != ii; i++)
+ {
+ std::cout<<" MuCCaFitter hits passed to computepam Z="
+ <<x[i]<<" phi="<<y[i]<<" zzz="<<z[i]<<" drift="<<r[i]<<" error="<<w[i]<<std::endl;
+ }
+ }
+ Fitter->Computeparam3(ii,z,y,r,sr);
+ if(m_debug){
+ std::cout<<" MuCCaFitter computed track a= "<<Fitter->get_a()
+ <<" da= "<<Fitter->get_da()
+ <<" b= "<<Fitter->get_b()
+ <<" db= "<<Fitter->get_db()
+ <<" corrab= "<<Fitter->get_corrab()
+ <<" chi2f= "<<Fitter->get_chi2f()
+ <<std::endl;
+ }
+
+ double afit=Fitter->get_a();
+ //double dafit=Fitter->get_da();
+ //double bfit=Fitter->get_b();
+ //double dbfit=Fitter->get_db();
+ double chi2f=Fitter->get_chi2f();
+
+ std::vector<double> dist(N);
+ std::vector<double> ddist(N);
+ double R,dis;
+ double theta = atan(afit);
+ if(theta<0.) theta= M_PI+theta;
+ double sinus = sin(theta);
+ double cosin = cos(theta);
+ double d = -( getZ( pos )-Zc)*sinus+( getY( pos )-Yc)*cosin;
+ if(m_debug){
+ std::cout << "MuCCaFitter>>> theta= " << theta << " sinus= " << sinus << " cosin= " << cosin << std::endl;
+ std::cout << "MuCCaFitter>>> getZ( pos )= " << getZ( pos ) << " getY( pos )= " << getY( pos ) << " Zc= " << Zc << " Yc= " << Yc<< std::endl;
+ }
+ double Szz(0),Syy(0),Szy(0),Syyzz(0),Att(0);
+ R=0;
+ for(int i=0;i<N;++i){
+ if(selection[i]) continue;
+ Szz += (z[i]-Zc)*(z[i]-Zc)*w[i];
+ Syy += (y[i]-Yc)*(y[i]-Yc)*w[i];
+ Szy += (y[i]-Yc)*(z[i]-Zc)*w[i];
+ Syyzz += ((y[i]-Yc)-(z[i]-Zc))*((y[i]-Yc)+(z[i]-Zc))*w[i];
+ dis = (y[i]-Yc)*cosin-(z[i]-Zc)*sinus;
+ if(dis>d){
+ R -= rw[i];
+ }else{
+ R += rw[i];
+ }
+ // std::cout << "MuCCaFitter>>> d= " << d <<" z[i]= "<<z[i] <<" y[i]= "<<y[i]<< " dis= " << dis << " R= " << R << " rw[i]= " << rw[i]<< std::endl;
+ }
+ Att = Syy + cosin*(2*sinus*Szy - cosin*Syyzz);
+ d = R/S;
+ if(m_debug)
+ std::cout << "MuCCaFitter>>> d= " << d << " R= "
+ << R << " S= " << S << " Att= " << Att << std::endl;
+ for(int i=0;i<N;++i){
+ if(selection[i]) continue;
+ dist[i] = cosin*(y[i]-Yc) - sinus*(z[i]-Zc) - d;
+ double dth = -(sinus*(y[i]-Yc) + cosin*(z[i]-Zc))*(sqrt(1./Att));
+ ddist[i] = sqrt( dth*dth + (1./S) );
+ }
+ if(m_debug) std::cout << "Transforming back to real world" << std::endl;
+ Amg::Vector3D ndir = getVec( 0., sinus, cosin );
+ Amg::Vector3D npos = getVec( 0., Yc + cosin*d, Zc - sinus*d );
+
+ if(m_debug) std::cout << "New line: position " << npos << " direction " << ndir << " chi2f " << chi2f << std::endl;
+
+ seg.set( chi2f/(N-2), npos, ndir );
+
+ MuonCalibSegment::MdtHitIt it = seg.mdtHOTBegin();
+
+ int i(0);
+ while(it!=seg.mdtHOTEnd()){
+ MdtCalibHitBase& hit = const_cast< MdtCalibHitBase& >( **it );
+
+ hit.setDistanceToTrack( dist[i], ddist[i] );
+ // std::cout << "MuCCaFitter>>> dist[i] " << dist[i] <<" dis= "<<dis<<" ddist= "<< ddist[i] << std::endl;
+ ++it;++i;
+ }
+
+ delete Fitter;
+ if(m_debug) std::cout << "fit done" << std::endl;
+ return true;
+}
+void MuCCaFitterImplementation::Computeparam3(int number_of_hits,std::vector<double> x,std::vector<double> y,std::vector<double> r,std::vector<double> sr)
+{
+/***************************************/
+/* Fit a line to n=number_of_hits hits */
+/***************************************/
+ double xout[100],yout[100];
+ // double grada[4],gradb[4];
+ double dist[100],rsub[100];
+ double chi2outn[4],chi2min;
+ // double chip[4];
+ double aref[4],bref[4];
+ // double sigmaa[4],sigmab[4];
+ int fhit,i,j,icouple=0,lhit;
+ // int segnobf[4];
+ // double da,db,hesseb,hessea,hesseab,chi2fin;
+ // for (int i = 0; i != number_of_hits; i++)
+ // {
+ // std::cout<<" MuCCaFitterImplementation hits passed to computepam Z="
+ // <<x[i]<<" phi="<<y[i]<<" drift="<<r[i]<<" error="<<sr[i]<<std::endl;
+ // }
+/* compute 4 tanget lines from first and last hit */
+ fhit = 0;
+ lhit = number_of_hits-1;
+ // Computelinparnew(x[fhit],y[fhit],r[fhit],sr[fhit],x[lhit],y[lhit],r[lhit],sr[lhit]);
+ Computelinparnew(x[fhit],y[fhit],r[fhit],x[lhit],y[lhit],r[lhit]);
+ for(int ii=0;ii<4;ii++){
+ bref[ii]=m_bfpar[ii];
+ aref[ii]=m_angularcoefficient[ii];
+ }
+/* choose the REFERENCE LINE (aref[icouple],bref[icouple])*/
+ for(i=0;i<4;i++){
+ for(j=0;j<number_of_hits;j++){
+ double d;
+ d = fabs(aref[i]*x[j]+bref[i]-y[j]);
+ dist[j] = d/sqrt(1.+aref[i]*aref[i]);
+ rsub[j] = r[j]-dist[j];
+ }
+ double chi2out = 0;
+ for(int ii=1;ii<number_of_hits-1;ii++){
+ chi2out += rsub[ii]*rsub[ii]/(sr[ii]*sr[ii]);
+ }
+ chi2outn[i]=chi2out;
+ }
+ chi2min = 999999999.;
+ for(i=0;i<4;i++){
+ if(chi2outn[i]<chi2min){
+ chi2min = chi2outn[i];
+ icouple = i;
+ }
+ }
+/* Compute TRACK POINTS */
+ Computehitpsemes(number_of_hits,x,y,r,aref[icouple],bref[icouple]);
+ for(i=0;i<number_of_hits;i++){
+ xout[i]=m_xpoint[i];
+ yout[i]=m_ypoint[i];
+ }
+/* Compute a & b parameters of the track from TRACK POINTS */
+ // int idim;
+ // ifail;
+ double aoutn, bout, sig2a, sig2b, corrab;
+ double temp,det;
+ double hesse[2][2];
+ double W,WX,WX2,WY,WY2,WXY;
+ // double errormatrix[2][2];
+ W=0.;
+ WX=0.;
+ WX2=0.;
+ WY=0.;
+ WY2=0.;
+ WXY=0;
+ for(int i=0;i<number_of_hits;i++){
+ temp = 1./(sr[i]*sr[i]);
+ W += temp;
+ WX += xout[i]*temp;
+ WX2+= xout[i]*xout[i]*temp;
+ WY += yout[i]*temp;
+ WY2+= yout[i]*yout[i]*temp;
+ WXY+= xout[i]*yout[i]*temp;
+ }
+ det = W*WX2-WX*WX;
+ aoutn = (W*WXY-WY*WX)/det;
+ bout = (WY*WX2-WX*WXY)/det;
+ hesse[1][1] = W;
+ hesse[0][0] = WX2;
+ hesse[1][0] = WX;
+ hesse[0][1] = WX;
+ // idim = 2;
+/* invert hessian matrix */
+ hesse[1][1] = hesse[0][0]/det;
+ hesse[0][0] = hesse[1][1]/det;
+ hesse[1][0] = -1.*(hesse[0][1]/det);
+ hesse[0][1] = -1.*(hesse[0][1]/det);
+ sig2a = hesse[0][0];
+ sig2b = hesse[1][1];
+ corrab = hesse[1][0];
+ double deno,btest,bs1,bs2,db1,db2;
+ btest = bout;
+ bout = 0;
+ deno = 0;
+ for(int i=0;i<number_of_hits;i++){
+ bs1 = (y[i]-aoutn*x[i]-r[i]*sqrt(1+aoutn*aoutn));
+ bs2 = (y[i]-aoutn*x[i]+r[i]*sqrt(1+aoutn*aoutn));
+ db1 = fabs(bs1-btest);
+ db2 = fabs(bs2-btest);
+ if(db1<db2){
+ bout += bs1/(sr[i]*sr[i]);
+ }else{
+ bout +=bs2/(sr[i]*sr[i]);
+ }
+ deno+= 1/(sr[i]*sr[i]);
+ }
+ bout /=deno;
+/* compute chi2 from residuals*/
+ double resd;
+ double chi2f = 0;
+ for(int i=0;i<number_of_hits;i++){
+ double xi,yi;
+ xi = (aoutn*y[i]-aoutn*bout+x[i])/(1.+aoutn*aoutn);
+ yi = aoutn*xi+bout;
+ resd = (sqrt((xi-x[i])*(xi-x[i])+(yi-y[i])*(yi-y[i]))-r[i])/sr[i];
+ chi2f += resd*resd;
+ }
+/* set global variables (method output)*/
+ m_aoutn=aoutn;
+ m_sig2a=sig2a;
+ m_bout=bout;
+ m_sig2b=sig2b;
+ m_corrab=corrab;
+ // if(fabs(atan(aoutn))>1.5) chi2f=100000000;
+ m_chi2f=chi2f;
+ // std::cout<<" MuCCaFitter::Computeparam3 computed track a= "<<m_aoutn
+ // <<" da= "<<m_sig2a
+ // <<" b= "<<m_bout
+ // <<" db= "<<m_sig2b
+ // <<" corrab= "<<m_corrab
+ // <<" chi2f= "<<m_chi2f
+ // <<std::endl;
+}
+
+double MuCCaFitterImplementation::get_a()
+{
+ return m_aoutn;
+}
+double MuCCaFitterImplementation::get_b() {return m_bout;}
+double MuCCaFitterImplementation::get_da() {return m_sig2a;}
+double MuCCaFitterImplementation::get_db() {return m_sig2b;}
+double MuCCaFitterImplementation::get_corrab() {return m_corrab;}
+double MuCCaFitterImplementation::get_chi2f() {return m_chi2f;}
+
+void MuCCaFitterImplementation::Computelinparnew(double x1, double y1, double r1, double x2, double y2, double r2)
+{
+/********************************************/
+/* Compute the 4 lines tangent to 2 circles */
+/********************************************/
+ double delta;
+ double averagephi,dist,dphiin,dphiex;
+ double bpar[4],phi;
+ double bfparn[2];
+ // int segnobfn[2];
+ // int segnoc[2][4],ncandid[4],ncand;
+ int ncandid[4],ncand;
+ double bcand[2][4];
+ int segnob[]={1,-1,1,-1};
+ int firsttime,i;
+ double angularcoefficient[4];
+ double bfpar[4];
+ // int segnobf[4];
+ double dy, dx;
+
+/* compute a parameters */
+ dx = x2-x1;
+ dy = y2-y1;
+ averagephi = atan2(dy,dx);
+ dist = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
+ delta = r2+r1;
+
+ dphiin = asin(delta/dist);
+
+ delta = r2-r1;
+ dphiex = asin(delta/dist);
+
+ int f=1;
+ phi = averagephi+f*dphiex;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[0] = tan(phi);
+
+ f=-1;
+ phi = averagephi+f*dphiex;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[1] = tan(phi);
+
+ f=1;
+ phi = averagephi+f*dphiin;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[2] = tan(phi);
+
+ f=-1;
+ phi = averagephi+f*dphiin;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[3] = tan(phi);
+
+/* Compute b parameters */
+ for(i=0;i<4;i++){
+ bpar[0] = y1-angularcoefficient[i]*x1+segnob[0]*r1*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ bpar[1] = y1-angularcoefficient[i]*x1+segnob[1]*r1*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ bpar[2] = y2-angularcoefficient[i]*x2+segnob[2]*r2*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ bpar[3] = y2-angularcoefficient[i]*x2+segnob[3]*r2*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ double delta = 0.00001;
+ ncand = 0;
+ if(fabs(bpar[0]-bpar[2])<delta){
+ bfparn[ncand]= bpar[0];
+ // segnobfn[ncand] = segnob[0];
+ ncand = ncand+1;
+ }
+ if(fabs(bpar[0]-bpar[3])<delta){
+ bfparn[ncand] = bpar[0];
+ // segnobfn[ncand] = segnob[0];
+ ncand = ncand+1;
+ }
+ if(fabs(bpar[1]-bpar[2])<delta){
+ bfparn[ncand] = bpar[1];
+ // segnobfn[ncand] = segnob[1];
+ ncand = ncand+1;
+ }
+ if(fabs(bpar[1]-bpar[3])<delta){
+ bfparn[ncand] = bpar[1];
+ // segnobfn[ncand] = segnob[1];
+ ncand = ncand+1;
+ }
+ bcand[0][i] = bfparn[0];
+ bcand[1][i] = bfparn[1];
+ // segnoc[0][i] = segnobfn[0];
+ // segnoc[1][i] = segnobfn[1];
+ ncandid[i] = ncand;
+ }
+ firsttime = 0;
+ for(i=0;i<4;i++){
+ if(ncandid[i]==1){
+ bfpar[i] = bcand[0][i];
+ // segnobf[i] = segnoc[0][i];
+ }else{
+ bfpar[i] = bcand[firsttime][i];
+ // segnobf[i] = segnoc[firsttime][i];
+ firsttime++;
+ }
+ }
+/* set global variables (method output)*/
+ for(i=0;i<4;i++){
+ m_bfpar[i]=bfpar[i];
+ m_angularcoefficient[i]=angularcoefficient[i];
+ }
+}
+
+void MuCCaFitterImplementation::Computelin(double x1, double y1, double r1, double x2, double y2, double r2)
+{
+/********************************************/
+/* Compute the 4 lines tangent to 2 circles */
+/********************************************/
+ double delta;
+ double averagephi,dist,dphiin,dphiex;
+ double bpar[4],phi;
+ // double distfp[4];
+ double bfparn[2];
+ double bcand[2][4];
+ // int segnobf[4],segnobfn[2];
+ // double sol[4];
+ int segnob[]={1,-1,1,-1};
+ // int segnoc[2][4],ncandid[4],ncand;
+ int ncandid[4],ncand;
+ int i,firsttime;
+ double angularcoefficient[4];
+ double bfpar[4];
+ double dy, dx;
+
+/* compute a parameters */
+ dx = x2-x1;
+ dy = y2-y1;
+ averagephi = atan2(dy,dx);
+ dist = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
+
+ delta = r2+r1;
+ dphiin = asin(delta/dist);
+
+ delta = r2-r1;
+ dphiex = asin(delta/dist);
+
+ int f=1;
+ phi = averagephi+f*dphiex;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[0] = tan(phi);
+
+ f=-1;
+ phi = averagephi+f*dphiex;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[1] = tan(phi);
+
+ f=1;
+ phi = averagephi+f*dphiin;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[2] = tan(phi);
+
+ f=-1;
+ phi = averagephi+f*dphiin;
+ if(phi < 0){phi = 6.2831853+(phi);}
+ angularcoefficient[3] = tan(phi);
+
+/* Compute b parameters */
+ for(i=0;i<4;i++){
+ bpar[0] = y1-angularcoefficient[i]*x1+segnob[0]*r1*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ bpar[1] = y1-angularcoefficient[i]*x1+segnob[1]*r1*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ bpar[2] = y2-angularcoefficient[i]*x2+segnob[2]*r2*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ bpar[3] = y2-angularcoefficient[i]*x2+segnob[3]*r2*
+ sqrt(1+angularcoefficient[i]*angularcoefficient[i]);
+ double delta = 0.00001;
+ ncand = 0;
+ if(fabs(bpar[0]-bpar[2])<delta){
+ bfparn[ncand] = bpar[0];
+ // segnobf[ncand] = segnob[0];
+ ncand = ncand+1;
+ }
+ if(fabs(bpar[0]-bpar[3])<delta){
+ bfparn[ncand] = bpar[0];
+ // segnobf[ncand] = segnob[0];
+ ncand = ncand+1;
+ }
+ if(fabs(bpar[1]-bpar[2])<delta){
+ bfparn[ncand] = bpar[1];
+ // segnobf[ncand] = segnob[1];
+ ncand = ncand+1;
+ }
+ if(fabs(bpar[1]-bpar[3])<delta){
+ bfparn[ncand] = bpar[1];
+ // segnobf[ncand] = segnob[1];
+ ncand = ncand+1;
+ }
+ bcand[0][i] = bfparn[0];
+ bcand[1][i] = bfparn[1];
+ // segnoc[0][i] = segnobfn[0];
+ // segnoc[1][i] = segnobfn[1];
+ ncandid[i] = ncand;
+ }
+ firsttime = 0;
+ for(i=0;i<4;i++){
+ if(ncandid[i]==1){
+ bfpar[i] = bcand[0][i];
+ // segnobf[i] = segnoc[0][i];
+ }else{
+ bfpar[i] = bcand[firsttime][i];
+ // segnobf[i] = segnoc[firsttime][i];
+ firsttime++;
+ }
+ }
+/* set global variables (method output)*/
+ for(int i=0; i<4; i++){
+ m_bfpar[i]=bfpar[i];
+ m_angularcoefficient[i]=angularcoefficient[i];
+ }
+}
+void MuCCaFitterImplementation::Computehitpsemes(int nhit,std::vector<double> xcirc, std::vector<double> ycirc, std::vector<double>rcirc, double a, double b)
+{
+/************************/
+/* Compute TRACK POINTS */
+/************************/
+ int i,nhitc;
+ // double xpoint[nhit],ypoint[nhit];
+/* loop over hit couples */
+ nhitc = nhit/2;
+ for(i=0;i<nhitc;i++){
+/* Compute TRACK POINTS for 2 circles */
+ computehitsfromcircles(xcirc[2*i],ycirc[2*i],rcirc[2*i],xcirc[2*i+1],ycirc[2*i+1],rcirc[2*i+1],a,b);
+ m_xpoint[2*i]=m_xout0;
+ m_ypoint[2*i]=m_yout0;
+ m_xpoint[2*i+1]=m_xout1;
+ m_ypoint[2*i+1]=m_yout1;
+ }
+ if(2*nhitc!=nhit){
+/* Compute TRACK POINTS for 2 circles */
+ computehitsfromcircles(xcirc[nhit-2],ycirc[nhit-2],rcirc[nhit-2],xcirc[nhit-1],ycirc[nhit-1],rcirc[nhit-1],a,b);
+ m_xpoint[nhit-2]=m_xout0;
+ m_ypoint[nhit-2]=m_yout0;
+ m_xpoint[nhit-1]=m_xout1;
+ m_ypoint[nhit-1]=m_yout1;
+ }
+}
+
+void MuCCaFitterImplementation::computehitsfromcircles(double x0,double y0,double r0,double x1,double y1,double r1,double a,double b)
+{
+/**************************************/
+/* Compute TRACK POINTS for 2 circles */
+/**************************************/
+ double ang[4],bpar[4];
+ int i;
+ double xout0=0,yout0=0;
+ double xout1=0,yout1=0;
+ double dist,dist1,xint,yint,xref,yref;
+
+/* Compute the 4 lines tangent to 2 circles */
+ Computelin(x0,y0,r0,x1,y1,r1);
+ for(i=0;i<4;i++){
+ bpar[i]=m_bfpar[i];
+ ang[i]=m_angularcoefficient[i];
+ }
+
+ dist = 9999999;
+/* for each of the 4 tangents : */
+ for(i=0;i<4;i++){
+/* compute tangent point */
+ xint = (x0+ang[i]*y0-ang[i]*bpar[i])/(1+ang[i]*ang[i]);
+ yint = ang[i]*(xint)+bpar[i];
+/* compute point from reference line */
+ double bprime;
+ if(a!=0){
+ bprime = y0+x0/a;
+ xref = (bprime-b)*a/(a*a+1);
+ yref = (b+bprime*a*a)/(a*a+1);
+ }else{
+ xref = x0;
+ yref = b;
+ }
+/* choose tangent point closest to that of the reference lines (TRACK POINT)*/
+ dist1 = sqrt((xint-xref)*(xint-xref)+(yint-yref)*(yint-yref));
+ if(dist1<dist){
+ dist = dist1;
+ xout0 = xint;
+ yout0 = yint;
+ }
+ }
+ dist = 9999999;
+ for(i=0;i<4;i++){
+/* compute tangent point */
+ xint = (x1+ang[i]*y1-ang[i]*bpar[i])/(1+ang[i]*ang[i]);
+ yint = ang[i]*(xint)+bpar[i];
+/* compute point from reference line */
+ double bprime;
+ if(a!=0){
+ bprime = y1+x1/a;
+ xref = (bprime-b)*a/(a*a+1);
+ yref = (b+bprime*a*a)/(a*a+1);
+ }else{
+ xref = x1;
+ yref = b;
+ }
+/* choose tangent point closest to that of the reference line (TRACK POINT)*/
+ dist1 = sqrt((xint-xref)*(xint-xref)+(yint-yref)*(yint-yref));
+ if(dist1<dist){
+ dist = dist1;
+ xout1 = xint;
+ yout1 = yint;
+ }
+ }
+/* set global variables (method output)*/
+ m_xout0=xout0;
+ m_yout0=yout0;
+ m_xout1=xout1;
+ m_yout1=yout1;
+}
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/QuasianalyticLineReconstruction.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/QuasianalyticLineReconstruction.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8de996c1286f3d73c7406a9fce58249c227b20e0
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/QuasianalyticLineReconstruction.cxx
@@ -0,0 +1,947 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 01.03.2006, AUTHOR: OLIVER KORTNER
+// Modified: 15.07.2006 by O. Kortner, error calculation corrected,
+// chi^2 refit functionality added.
+// 13.01.2007 by O. Kortner, bug fix in candidate treatment, some
+// candidates were considered to be
+// identical although different;
+// modifications to improve the
+// reconstruction efficiency at very high
+// background count rates.
+// 27.03.2007 by O. Kortner, distances with signs filled into
+// MuonCalibSegment.
+// 23.03.2007 by O. Kortner, isnan check for chi^2.
+// 08.06.2007 by O. Kortner, final track segment has rphi position and
+// direction of the initial segment.
+// 23.06.2006 by O. Kortner, add convention for rphi track if the
+// pattern recognition has failed to
+// provide it.
+// 26.11.2007 by O. Kortner, fix for segment refinement.
+// 13.12.2007 by O. Kortner, time-out added.
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibFitters/QuasianalyticLineReconstruction.h"
+#include <iostream>
+#include <fstream>
+#include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+#include "time.h"
+#include <cmath>
+
+//:::::::::::::::::::::::
+//:: NAMESPACE SETTING ::
+//:::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS ::
+//:: QuasianalyticLineReconstruction ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD m_init ::
+//:::::::::::::::::::
+
+void QuasianalyticLineReconstruction::m_init(void) {
+
+ m_init(0.5*CLHEP::mm); // default road width = 0.5 CLHEP::mm
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD m_init(const double & r_road_width) ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::
+
+void QuasianalyticLineReconstruction::m_init(const double & r_road_width) {
+
+//:::::::::::::::
+//:: VARIABLES ::
+//:::::::::::::::
+
+ Amg::Vector3D null_vec(0.0, 0.0, 0.0); // auxiliary 0 vector
+
+//::::::::::::::::::
+//:: SET TIME-OUT ::
+//::::::::::::::::::
+
+ m_time_out = 10.0;
+
+//::::::::::::::::::::::::::::
+//:: SET THE MAXIMUM RADIUS ::
+//::::::::::::::::::::::::::::
+
+ m_r_max = 15.0;
+
+//::::::::::::::::::::::::
+//:: SET THE ROAD WIDTH ::
+//::::::::::::::::::::::::
+
+ m_road_width = r_road_width;
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: INITIALIZE PRIVATE VARIABLES WHICH ARE ACCESSIBLE BY METHODS ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+ m_nb_track_hits = 0;
+ m_chi2 = 0.0;
+ m_track = MTStraightLine(null_vec, null_vec, null_vec, null_vec);
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: SET THE TRACK IN THE x1-x3 PLANE (LATER VERSION MAY ALLOW TO SET ::
+//:: USER-DEFINED VALUES) ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+ m_m_x1 = 1.0; m_b_x1 = 0.0;
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::
+//:: METHOD tangent ::
+//::::::::::::::::::::
+
+MTStraightLine QuasianalyticLineReconstruction::tangent(
+ const Amg::Vector3D & r_w1,
+ const double & r_r1, const double & r_sigma12,
+ const Amg::Vector3D & r_w2, const double & r_r2,
+ const double & r_sigma22, const int & r_case) const {
+
+//:::::::::::::::
+//:: VARIABLES ::
+//:::::::::::::::
+
+ double sinalpha; // auxiliary sinus of an angle
+ double cosalpha; // auxiliary sinus of an angle
+ Amg::Vector3D e2prime(0.,0.,0.), e3prime(0.,0.,0.); // auxiliary direction vectors
+ MTStraightLine tang; // tangent to drift circles of a hit pair
+ Amg::Vector3D p1(0.,0.,0.), p2(0.,0.,0.); // hit points defining a tangent
+ Amg::Vector3D null_vec(0.0, 0.0, 0.0); // auxiliary 0 vector
+ double mx1, bx1, mx2, bx2; // auxiliary track parameters
+
+//::::::::::::::::::::::::::::::::::::::::::::
+//:: CHECK WHETHER THE SELECTED CASE EXISTS ::
+//::::::::::::::::::::::::::::::::::::::::::::
+
+ if (r_case<1 || r_case>4) {
+ cerr << endl
+ << "Class QuasianalyticLineReconstruction, "
+ << "method tangent: ERROR!\n"
+ << "Illegal case " << r_case << ", must be 1,2,3, or 4."
+ << endl;
+ exit(1);
+ }
+
+//:::::::::::::::::::::::::::::::::::::
+//:: CALCULATE THE REQUESTED TANGENT ::
+//:::::::::::::::::::::::::::::::::::::
+
+// local coordinate axis vectors //
+ e3prime = (r_w2-r_w1).unit();
+ e2prime = Amg::Vector3D(0.0, e3prime.z(), -e3prime.y());
+
+// case 1 and 2 //
+ if (r_case==1 || r_case==2) {
+
+ sinalpha = fabs(r_r2-r_r1)/(r_w2-r_w1).mag();
+ cosalpha = sqrt(1.0-sinalpha*sinalpha);
+
+ // case 1 //
+ if (r_case==1) {
+ p1 = r_w1+((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r1*
+ (cosalpha*e2prime-sinalpha*e3prime);
+ p2 = r_w2+((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r2*
+ (cosalpha*e2prime-sinalpha*e3prime);
+ }
+
+ // case 2 //
+ if (r_case==2) {
+ p1 = r_w1-((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r1*
+ (cosalpha*e2prime+sinalpha*e3prime);
+ p2 = r_w2-((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r2*
+ (cosalpha*e2prime+sinalpha*e3prime);
+ }
+
+
+ }
+
+// case 3 and 4 //
+ if (r_case==3 || r_case==4) {
+
+ sinalpha = (r_r1+r_r2)/(r_w2-r_w1).mag();
+ cosalpha = sqrt(1.0-sinalpha*sinalpha);
+
+ // case 3 //
+ if (r_case==3) {
+// p1 = r_w1+r_r1*(cosalpha*e2prime-sinalpha*e3prime);
+// p2 = r_w2-r_r2*(cosalpha*e2prime-sinalpha*e3prime);
+ p1 = r_w1+r_r1*(cosalpha*e2prime+sinalpha*e3prime);
+ p2 = r_w2-r_r2*(cosalpha*e2prime+sinalpha*e3prime);
+ }
+
+ // case 4 //
+ if (r_case==4) {
+ p1 = r_w1-r_r1*(cosalpha*e2prime-sinalpha*e3prime);
+ p2 = r_w2+r_r2*(cosalpha*e2prime-sinalpha*e3prime);
+ }
+
+ }
+
+// calculation of the tangent and estimation of its errors //
+ if ((p2-p1).z()!=0.0) {
+ tang = MTStraightLine(p1, p2-p1, null_vec, null_vec);
+ } else {
+ Amg::Vector3D direction(p2-p1);
+ direction[2] = 1.0e-99;
+ tang = MTStraightLine(p1, direction, null_vec, null_vec);
+ }
+ mx1 = tang.m_x1(); bx2 = tang.b_x1();
+ mx2 = tang.m_x2(); bx2 = tang.b_x2();
+ tang = MTStraightLine(mx1, bx1, mx2, bx2,1.0, 1.0,
+ sqrt(r_sigma12+r_sigma22)/fabs(p2.z()-p1.z()),
+ sqrt(r_sigma12+p1.z()*p1.z()*(r_sigma12+r_sigma22)/
+ std::pow(p2.z()-p1.z(), 2)));
+ // errors in mx1 and bx1 are arbitrary since they are
+ // not used at a later stage.
+
+
+ return tang;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::
+//:: METHOD m_track_candidate ::
+//::::::::::::::::::::::::::::::
+
+MTStraightLine QuasianalyticLineReconstruction::m_track_candidate(
+ const IndexSet & r_index_set,
+ const int & r_k_cand,
+ const int & r_l_cand,
+ const int & r_cand_case,
+ vector<Amg::Vector3D> r_w,
+ vector<double> r_r,
+ vector<double> r_sigma2,
+ double & r_chi2) const {
+
+//:::::::::::::::
+//:: VARIABLES ::
+//:::::::::::::::
+
+ int nb_hits(r_index_set.size()); // number of hits used in the track
+ // reconstruction
+ int nb_tangents(0); // number of tangents used in the track
+ // reconstruction
+ vector<double> mx1, bx1, mx2, bx2; // slopes and intercepts of the
+ // tangents building the track
+ vector<double> mx1_err, bx1_err, mx2_err, bx2_err; // their errors
+ MTStraightLine aux_track; // auxiliary straight line
+ MTStraightLine tang; // tangent to drift circles of a hit pair
+ Amg::Vector3D d1(0.0, 0.0, 0.0), d2(0.0, 0.0, 0.0); // auxiliary distance vectors
+ Amg::Vector3D e2prime(0.0, 0.0, 0.0), e3prime(0.0, 0.0, 0.0); // auxiliary basis vectors
+ Amg::Vector3D a(0.0, 0.0, 0.0), u(0.0, 0.0, 0.0); // position and direction vector of the candidate
+ // tangent
+ double sinalpha; // auxiliary sinus of an angle
+ double cosalpha; // auxiliary sinus of an angle
+ Amg::Vector3D p1(0.0, 0.0, 0.0), p2(0.0, 0.0, 0.0); // hit points defining a tangent
+ Amg::Vector3D null_vec(0.0, 0.0, 0.0); // auxiliary 0 vector
+ double sum[4]; // auxiliary summation variables
+
+//:::::::::::::::::::::
+//:: GET THE TANGENT ::
+//:::::::::::::::::::::
+
+// calculate the tangent //
+ aux_track = tangent(r_w[r_k_cand], r_r[r_k_cand], r_sigma2[r_k_cand],
+ r_w[r_l_cand], r_r[r_l_cand], r_sigma2[r_l_cand],
+ r_cand_case);
+ a = aux_track.positionVector();
+ u = (aux_track.directionVector()).unit();
+
+// store the parameters of this tangent //
+ mx1.push_back(aux_track.m_x1());
+ mx1_err.push_back(aux_track.m_x1_error());
+ bx1.push_back(aux_track.b_x1());
+ bx1_err.push_back(aux_track.b_x1_error());
+ mx2.push_back(aux_track.m_x2());
+ mx2_err.push_back(aux_track.m_x2_error());
+ bx2.push_back(aux_track.b_x2());
+ bx2_err.push_back(aux_track.b_x2_error());
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: LOOP OVER THE OTHER HITS AND CALCULATE TANGENTS :: //
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+ for (int kk=0; kk<nb_hits-1; kk++) {
+ for (int ll=kk+1; ll<nb_hits; ll++) {
+
+ int k(r_index_set[kk]);
+ int l(r_index_set[ll]);
+
+ if (k==r_k_cand && l==r_l_cand) {
+ continue;
+ }
+
+// local coordinate axis vectors //
+ e3prime = (r_w[l]-r_w[k]).unit();
+ e2prime = Amg::Vector3D(0.0, e3prime.z(), -e3prime.y());
+
+// distance vectors //
+ d1 = (r_w[k]-a).dot(u)*u-(r_w[k]-a); d1[0] = (0.0);
+ d2 = (r_w[l]-a).dot(u)*u-(r_w[l]-a); d2[0] = (0.0);
+
+// one must distinguish 2 cases //
+
+ // case 1: candidate tangent passes both wires on the same side //
+ if (d1.dot(d2)>=0) {
+
+ sinalpha = fabs(r_r[l]-r_r[k])/(r_w[l]-r_w[k]).mag();
+ cosalpha = sqrt(1.0-sinalpha*sinalpha);
+
+ if (d1.dot(e2prime)>=0) {
+ if (r_r[k]<=r_r[l]) {
+ p1 = r_w[k]+r_r[k]*(cosalpha*e2prime-
+ sinalpha*e3prime);
+ } else {
+ p1 = r_w[k]+r_r[k]*(cosalpha*e2prime+
+ sinalpha*e3prime);
+ }
+ }
+ else {
+ if (r_r[k]>=r_r[l]) {
+ p1 = r_w[k]-r_r[k]*(cosalpha*e2prime-
+ sinalpha*e3prime);
+ } else {
+ p1 = r_w[k]-r_r[k]*(cosalpha*e2prime+
+ sinalpha*e3prime);
+ }
+ }
+
+ if (d2.dot(e2prime)>=0) {
+ if (r_r[k]<=r_r[l]) {
+ p2 = r_w[l]+r_r[l]*(cosalpha*e2prime-
+ sinalpha*e3prime);
+ } else {
+ p2 = r_w[l]+r_r[l]*(cosalpha*e2prime+
+ sinalpha*e3prime);
+ }
+ }
+ else {
+ if (r_r[k]>=r_r[l]) {
+ p2 = r_w[l]-r_r[l]*(cosalpha*e2prime-
+ sinalpha*e3prime);
+ } else {
+ p2 = r_w[l]-r_r[l]*(cosalpha*e2prime+
+ sinalpha*e3prime);
+ }
+ }
+
+ }
+
+ // case 2: candidate tangent passes both wires on opposite sides //
+ if (d1.dot(d2)<0) {
+
+ sinalpha = (r_r[l]+r_r[k])/(r_w[l]-r_w[k]).mag();
+ cosalpha = sqrt(1.0-sinalpha*sinalpha);
+
+ // case 2(a) //
+ if (d1.dot(e2prime)>=0 && d2.dot(e2prime)<=0) {
+ p1 = r_w[k]+r_r[k]*(cosalpha*e2prime+
+ sinalpha*e3prime);
+ p2 = r_w[l]-r_r[l]*(cosalpha*e2prime+
+ sinalpha*e3prime);
+ }
+
+ // case 2(b) //
+ if (d1.dot(e2prime)<0 && d2.dot(e2prime)>=0) {
+ p1 = r_w[k]-r_r[k]*(cosalpha*e2prime-
+ sinalpha*e3prime);
+ p2 = r_w[l]+r_r[l]*(cosalpha*e2prime-
+ sinalpha*e3prime);
+ }
+
+ }
+
+// get the slope and intercept of the new tangents and error estimates //
+ if ((p2-p1).z()!=0.0) {
+ tang = MTStraightLine(p1, p2-p1, null_vec, null_vec);
+ } else {
+ Amg::Vector3D direction(p2-p1);
+ direction[2] = (1.0e-99);
+ tang = MTStraightLine(p1, direction,
+ null_vec, null_vec);
+ }
+ mx1.push_back(tang.m_x1());
+ bx1.push_back(tang.b_x1());
+ mx2.push_back(tang.m_x2());
+ bx2.push_back(tang.b_x2());
+
+ mx1_err.push_back(1.0);
+ bx1_err.push_back(1.0);
+
+ mx2_err.push_back((r_sigma2[k]+r_sigma2[l])/
+ std::pow(p2.z()-p1.z(), 2));
+ bx2_err.push_back(r_sigma2[k]+
+ p1.z()*p1.z()*(r_sigma2[k]+r_sigma2[l])/
+ std::pow(p2.z()-p1.z(),2));
+ // errors in x1 and x2 are arbitrary since they are not
+ // used at a later stage.
+
+// increase the number of tangents //
+ nb_tangents = nb_tangents+1;
+
+ }
+ }
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: CALCULATE A WEIGHTED AVERAGE OVER THE TANGENTS AND THE chi^2 FOR ::
+//:: THE RECONSTRUCTED TRAJECTORY ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+ sum[0] = 0.0; sum[1] = 0.0; sum[2] = 0.0; sum[3] = 0.0;
+ for (int k=0; k<nb_tangents; k++) {
+
+ if (mx2_err[k]>0) {
+ sum[0] = sum[0] + mx2[k]/mx2_err[k];
+ sum[1] = sum[1] + 1.0/mx2_err[k];
+ } else {
+ sum[0] = sum[0] + mx2[k];
+ sum[1] = sum[1] + 1.0;
+ }
+ if (bx2_err[k]>0) {
+ sum[2] = sum[2] + bx2[k]/bx2_err[k];
+ sum[3] = sum[3] + 1.0/bx2_err[k];
+ } else {
+ sum[2] = sum[2] + bx2[k];
+ sum[3] = sum[3] + 1.0;
+ }
+
+ }
+
+ aux_track = MTStraightLine(0.0, 0.0, sum[0]/sum[1], sum[2]/sum[3],
+ 0.0, 1.0, sqrt(1.0/sum[1]), sqrt(1.0/sum[3]));
+
+ r_chi2 = 0.0;
+ for (int k=0; k<nb_hits; k++) {
+ r_chi2 = r_chi2+std::pow(aux_track.distFromLine(
+ r_w[r_index_set[k]])
+ -r_r[r_index_set[k]], 2)/
+ r_sigma2[k];
+ }
+
+ return aux_track;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD roadWidth ::
+//::::::::::::::::::::::
+
+double QuasianalyticLineReconstruction::roadWidth(void) const {
+
+ return m_road_width;
+
+}
+
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::
+//:: METHOD numberOfTrackHits ::
+//::::::::::::::::::::::::::::::
+
+unsigned int QuasianalyticLineReconstruction::numberOfTrackHits(void) const {
+
+ return m_nb_track_hits;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD trackHits ::
+//::::::::::::::::::::::
+
+const vector<const MdtCalibHitBase*> & QuasianalyticLineReconstruction::trackHits(void) const {
+
+ return m_track_hits;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD chi2 ::
+//:::::::::::::::::
+
+double QuasianalyticLineReconstruction::chi2(void) const {
+
+ if (m_nb_track_hits == 0) {
+ return -1;
+ }
+ return m_chi2;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::
+//:: METHOD chi2PerDegreesOfFreedom ::
+//::::::::::::::::::::::::::::::::::::
+
+double QuasianalyticLineReconstruction::chi2PerDegreesOfFreedom(
+ void) const{
+
+ if (m_nb_track_hits<3) {
+ return -1;
+ }
+ return chi2()/static_cast<double>(numberOfTrackHits()-2);
+
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::
+//:: METHOD track ::
+//::::::::::::::::::
+
+MTStraightLine QuasianalyticLineReconstruction::track(void) const {
+
+ return m_track;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD setRoadWidth ::
+//:::::::::::::::::::::::::
+
+void QuasianalyticLineReconstruction::setRoadWidth(
+ const double & r_road_width) {
+
+ m_road_width = fabs(r_road_width);
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD setTimeOut ::
+//:::::::::::::::::::::::
+
+void QuasianalyticLineReconstruction::setTimeOut(const double & time_out) {
+
+ m_time_out = time_out;
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD setTimeOut ::
+//:::::::::::::::::::::::
+
+void setTimeOut(const double & /*time_out*/) {
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD fit(.) ::
+//:::::::::::::::::::
+
+bool QuasianalyticLineReconstruction::fit(MuonCalibSegment & r_segment) const {
+
+// select all hits //
+ HitSelection selection(r_segment.mdtHitsOnTrack(), 0);
+
+// call the other fit function //
+ return fit(r_segment, selection);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD fit(.,.) ::
+//:::::::::::::::::::::
+
+bool QuasianalyticLineReconstruction::fit(MuonCalibSegment & r_segment,
+ HitSelection r_selection) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ time_t start,end; // start and end time (needed for time-out)
+ time(&start);
+ double diff; // difference of start and end time (needed for time-out)
+ unsigned int nb_selected_hits(0); // number of selected hits
+ vector<const MdtCalibHitBase*> selected_hits; // vector of pointers to the
+ // selected hits
+ vector<unsigned> selected_hits_index(r_selection.size());
+ // vector containing the indices
+ // of the selected hits (needed
+ // a requested refit)
+ vector<Amg::Vector3D> w; // wire coordinates
+ vector<double> r; // drift CLHEP::radii of the selected hits
+ vector<double> sigma2; // sigma(r)^2 (spatial resolution in r)
+ int counter1, counter2, counter3; // auxiliary counters
+ int max_cand_hits; // the maximum number of hits on a track candidate
+ // found so far
+ int max_cand_HOTs; // the maximum number of hit tubes on a track
+ // candidate found so far
+ int nb_candidates; // number of candidate tracks
+ vector<int> k_cand, l_cand; // candidate track index
+ vector<int> cand_case; // tangent case for the candidate
+ vector<int> nb_HOTs; // number of hits on a candidate
+ int candidate; // index of the candidate which becomes the track
+ IndexSet aux_set; // auxiliary index set
+ vector<IndexSet> index_set; // index set for the storage of the hits
+ // forming a track
+ vector<double> intercept; // intercepts of track candidates
+ MTStraightLine tangents[4]; // the four tangents to the drift circles of a
+ // hit pair
+ MTStraightLine aux_track; // auxiliary track
+ double aux_chi2; // auxiliary chi^2 variable
+ Amg::Vector3D null(0.0, 0.0, 0.0);
+ Amg::Vector3D xhat(1.0, 0.0, 0.0);
+ MTStraightLine initial_track(r_segment.position(),
+ r_segment.direction(), null, null);
+ // initial track stored in the segment
+ Amg::Vector3D aux_pos, aux_dir; // auxiliary position and direction vectors
+
+/////////////////////////////////////
+// GET THE NUMBER OF SELECTED HITS //
+/////////////////////////////////////
+ if (r_segment.mdtHitsOnTrack() != r_selection.size()) {
+ cerr << "\n"
+ << "Class QuasianalyticLineReconstruction, "
+ << "METHOD fit(., .): WARNING!\n"
+ << "Vector with selected hits unequal to the number "
+ << "of hits on the segment!\n"
+ << "The user selection will be ignored!\n";
+ r_selection.clear();
+ r_selection.assign(r_segment.hitsOnTrack(), 0);
+ } else {
+ for (unsigned int k=0; k<r_selection.size(); k++) {
+ if (r_selection[k] == 0) {
+ nb_selected_hits = nb_selected_hits+1;
+ }
+ }
+ }
+
+//////////////////////////////////////////////////////////////////////////
+// FIX POTENTIAL SECOND-COORDINATE PROBLEM OF THE INITIAL TRACK SEGMENT //
+//////////////////////////////////////////////////////////////////////////
+
+ if (r_segment.direction().z()==0.0) {
+ Amg::Vector3D dir(r_segment.direction().x(),
+ r_segment.direction().y(),
+ 1.0);
+ initial_track = MTStraightLine(r_segment.position(), dir,
+ null, null);
+ }
+
+///////////////////////////////////////////////////////////////////////////
+// RETURN, IF THERE ARE LESS THAN 3 HITS, I.E. THE TRACK IS NOT UNIQUELY //
+// DEFINED BY THE HITS. //
+///////////////////////////////////////////////////////////////////////////
+
+ if (nb_selected_hits<3) {
+ cerr << "\n"
+ << "Class QuasianalyticLineReconstruction, "
+ << "METHOD fit(., .): WARNING!\n"
+ << "Too few hits for the track reconstructions!\n";
+ return false;
+ }
+
+/////////////////////////////////////////////////
+// COPY THE WIRE COORDINATES TO A LOCAL VECTOR //
+/////////////////////////////////////////////////
+
+// vector assignments //
+ selected_hits.clear();// = vector<const MdtCalibHitBase*>(nb_selected_hits);
+ w.clear();// = vector<Amg::Vector3D>(nb_selected_hits);
+ r.clear();// = vector<double>(nb_selected_hits);
+ sigma2.clear();// = vector<double>(nb_selected_hits);
+ selected_hits_index.clear();//
+// vector filling //
+ counter2 = 0;
+ for (unsigned int k=0; k<r_segment.mdtHitsOnTrack(); k++) {
+
+ // skip the hit, if it has not been selected //
+ const MdtCalibHitBase *hit = (r_segment.mdtHOT())[k];
+ if (r_selection[k]==1 || hit->sigmaDriftRadius()>100.0) {
+ continue;
+ }
+
+ // copy the hit information into local vectors //
+ selected_hits.push_back(hit);//[counter2] = hit;
+ selected_hits_index.push_back(k);//[counter2] = k;
+ w.push_back(Amg::Vector3D(0.0, (hit->localPosition()).y(), (hit->localPosition()).z()));//[counter2] = Amg::Vector3D(0.0, (hit->localPosition()).y(), (hit->localPosition()).z());
+ r.push_back(fabs(hit->driftRadius()));//[counter2] = fabs(hit->driftRadius());
+ sigma2.push_back(hit->sigma2DriftRadius());//[counter2] = hit->sigma2DriftRadius();
+ // if the spatial resolution has not been set, set it to 0.1 CLHEP::mm //
+ if (sigma2[counter2] == 0) {
+ sigma2[counter2] = std::pow(0.1*CLHEP::mm, 2);
+ }
+
+ counter2 = counter2+1;
+
+ }
+ nb_selected_hits=selected_hits.size();
+ if(nb_selected_hits<3)
+ return false;
+//////////////////////////
+// TRACK RECONSTRUCTION //
+//////////////////////////
+
+// reset counters //
+ max_cand_hits = 0;
+ max_cand_HOTs = 0;
+ nb_candidates = 0;
+// 1st step: loop over all hit pairs and determine a track candidate //
+ for (unsigned int k=0; k<nb_selected_hits; k++) {
+ for (unsigned int l=k+1; l<nb_selected_hits; l++) {
+
+// time-out //
+ time (&end);
+ diff = difftime (end,start);
+ if (diff>m_time_out) {
+ cerr << endl
+ << "Class QuasianalyticLineReconstruction: "
+ << "time-out for track finding after "
+ << m_time_out
+ << " seconds!\n";
+ return false;
+ }
+
+// 2nd step: determine all four tangents to the drift circles of the hit pair //
+ for (unsigned int r_case=1; r_case<5; r_case++) {
+ tangents[r_case-1] = tangent(w[l], r[l], sigma2[l],
+ w[k], r[k], sigma2[k], r_case);
+ }
+
+// 3rd step: determine additional track points within the road width around //
+// the four tangents //
+ for (unsigned int r_case=1; r_case<5; r_case++) {
+ bool same(false);
+ counter1 = 0;
+ counter3 = 0;
+ vector<int> indices; // indices of the hits forming a
+ // track
+ for (unsigned int n=0; n<nb_selected_hits; n++) {
+ MTStraightLine aux_line(w[n], xhat, null, null);
+ if (fabs(fabs(tangents[r_case-1].signDistFrom(
+ aux_line))-r[n])<m_r_max) {
+ counter3++;
+ }
+ if (fabs(fabs(tangents[r_case-1].signDistFrom(
+ aux_line))-
+ r[n])<=m_road_width) {
+ counter1 = counter1+1;
+ indices.push_back(n);
+ }
+ }
+ if (counter1>2) {
+ aux_set = IndexSet(counter1, indices);
+ aux_set.sort();
+ for (int ca=0; ca<nb_candidates; ca++) {
+ if (aux_set==index_set[ca] &&
+ fabs(intercept[ca]-
+ tangents[r_case-1].b_x2())
+ <m_road_width) {
+ same = true;
+ break;
+ }
+ }
+ if (same) {
+ continue;
+ }
+ nb_candidates = nb_candidates+1;
+ k_cand.push_back(k);
+ l_cand.push_back(l);
+ cand_case.push_back(r_case);
+ index_set.push_back(aux_set);
+ intercept.push_back(tangents[r_case-1].b_x2());
+ nb_HOTs.push_back(counter3);
+ if (counter1>max_cand_hits) {
+ max_cand_hits = counter1;
+ }
+ }
+ }
+
+ }
+ }
+
+// 4th step: reconstruct a straight-line trajectory for all candidates with //
+// the maximum number of hits, keep the one with the smallest chi^2 //
+ m_nb_track_hits = max_cand_hits;
+ m_chi2 = -1.0;
+ if (m_nb_track_hits<3) {
+ return false;
+ }
+
+ candidate = 0;
+ for (int ca=0; ca<nb_candidates; ca++) {
+ if (index_set[ca].size()<static_cast<unsigned int>(
+ m_nb_track_hits)) {
+ continue;
+ }
+
+ aux_track = m_track_candidate(index_set[ca],
+ k_cand[ca], l_cand[ca],
+ cand_case[ca], w, r, sigma2, aux_chi2);
+ if (nb_HOTs[ca]>=max_cand_HOTs) {
+ max_cand_HOTs = nb_HOTs[ca];
+ if (m_chi2==-1.0 || m_chi2>aux_chi2) {
+ m_chi2 = aux_chi2;
+ m_track = aux_track;
+ candidate = ca;
+ }
+ }
+ }
+
+ m_track_hits.resize(m_nb_track_hits);
+ for (int k=0; k<m_nb_track_hits; k++) {
+ m_track_hits[k] = selected_hits[index_set[candidate][k]];
+ }
+
+ // make the segment in rphi as the initial segment //
+ aux_pos = Amg::Vector3D(initial_track.b_x1(), m_track.b_x2(), 0);
+ aux_dir = Amg::Vector3D(initial_track.m_x1(), m_track.m_x2(), 1.0);
+ m_track = MTStraightLine(aux_pos, aux_dir, null, null);
+
+// 5th step: update the segment//
+ if (std::isnan(m_chi2)) {
+ m_chi2=1.0e6;
+ }
+ r_segment.set(m_chi2/static_cast<double>(m_nb_track_hits-2),
+ m_track.positionVector(), m_track.directionVector());
+
+
+ vector<unsigned int> refit_hit_selection(r_selection.size(), 1);
+ for (int k=0; k<m_nb_track_hits; k++) {
+ refit_hit_selection[selected_hits_index[
+ index_set[candidate][k]]] = 0;
+ }
+ if(!m_refit && m_refine_segment) {
+ r_segment.refineMdtSelection(refit_hit_selection);
+ }
+
+ // chi^2 refit, if requested //
+ if (m_refit) {
+ if (m_nfitter.fit(r_segment, refit_hit_selection)) {
+// m_chi2 = r_segment.chi2()*(m_nb_track_hits-2);
+ Amg::Vector3D dir(r_segment.direction());
+ if (dir.z()==0.0) {
+ dir[2] = (1.0e-99);
+ }
+ MTStraightLine aux_line(r_segment.position(),
+ dir,
+ Amg::Vector3D(0.0, 0.0, 0.0),
+ Amg::Vector3D(0.0, 0.0, 0.0));
+ double m_err(m_track.m_x2_error());
+ double b_err(m_track.b_x2_error());
+ m_track = MTStraightLine(0.0, 0.0,
+ aux_line.m_x2(),
+ aux_line.b_x2(),
+ 0.0, 0.0,
+ m_err, b_err);
+ aux_pos = Amg::Vector3D(initial_track.b_x1(),
+ m_track.b_x2(), 0);
+ aux_dir = Amg::Vector3D(initial_track.m_x1(),
+ m_track.m_x2(), 1.0);
+ m_track = MTStraightLine(aux_pos, aux_dir, null, null);
+ // recompute chi^2 because of different convention in DCSLFitter //
+ m_chi2 = 0.0;
+ for (int k=0; k<m_nb_track_hits; k++) {
+ MTStraightLine wire(Amg::Vector3D(0.0,
+ m_track_hits[k]->localPosition().y(),
+ m_track_hits[k]->localPosition().z()),
+ xhat, null, null);
+ double d(fabs(m_track.signDistFrom(wire)));
+ m_chi2 = m_chi2+
+ std::pow(d-m_track_hits[k]->driftRadius(),
+ 2)/
+ m_track_hits[k]->sigma2DriftRadius();
+ }
+ r_segment.set(
+ m_chi2/static_cast<double>(m_nb_track_hits-2),
+ m_track.positionVector(), m_track.directionVector());
+ if (std::isnan(m_chi2)) {
+ m_chi2=1.0e6;
+ }
+ if (m_refine_segment) {
+ r_segment.refineMdtSelection(
+ refit_hit_selection);
+ }
+ return true;
+ }
+ return false;
+ }
+
+ // finish the update (necessary, if no refit has been performed) //
+ MuonCalibSegment::MdtHitIt it = r_segment.mdtHOTBegin();
+ while(it!=r_segment.mdtHOTEnd()){
+ MdtCalibHitBase& hit = const_cast< MdtCalibHitBase& >( **it );
+
+ Amg::Vector3D pos(0.0, (hit.localPosition()).y(),
+ (hit.localPosition()).z());
+ // wire position
+ MTStraightLine aux_line(pos, xhat, null, null);
+// double dist(fabs(m_track.signDistFrom(aux_line))); // track distance
+ double dist(m_track.signDistFrom(aux_line)); // track distance
+ double dist_err(sqrt(std::pow(pos.z()*m_track.m_x2_error(), 2)+
+ std::pow(m_track.b_x2_error(), 2)));
+ // approximate error of the track distance
+ hit.setDistanceToTrack(dist, dist_err);
+
+ ++it;
+ }
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD printLevel ::
+//:::::::::::::::::::::::
+
+void QuasianalyticLineReconstruction::printLevel(int level) {
+
+ cerr << "\n"
+ << "Class QuasianalyticLineReconstruction, "
+ << "method printLevel: WARNING!\n"
+ << "Print level " << level << " is ignored.\n";
+ return;
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/StraightPatRec.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/StraightPatRec.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..23c19e83a5b7ea2ba7b3ea4b054440df02d8cf2a
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibFitters/src/StraightPatRec.cxx
@@ -0,0 +1,797 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+// 27.10.2007, AUTHOR: OLIVER KORTNER
+// Modified: 26.11.2007 by O. Kortner, fix for segment refinement.
+// 13.12.2007 by O. Kortner, time-out added.
+// 07.08.2008 by O. Kortner, bug fix when hits are disabled.
+// 19.11.2008 by I. Potrap: 1) several bugs have been fixed;
+// 2) external fitter(DCLS) replaced by
+// internal code of linear-least-chi2-fit
+// 01.12.2008 by I. Potrap, track() method and segment position are fixed
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS ::
+//:: StraightPatRec ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+//::::::::::::::::::
+//:: HEADER FILES ::
+//::::::::::::::::::
+
+#include "MdtCalibFitters/StraightPatRec.h"
+#include <iostream>
+#include <fstream>
+#include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+
+#include "MuonCalibMath/Combination.h"
+#include "time.h"
+
+#include "EventPrimitives/EventPrimitives.h"
+#include "GeoPrimitives/GeoPrimitivesHelpers.h"
+
+//:::::::::::::::::::::::
+//:: NAMESPACE SETTING ::
+//:::::::::::::::::::::::
+
+using namespace MuonCalib;
+using namespace std;
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD m_init ::
+//:::::::::::::::::::
+
+void StraightPatRec::m_init(void) {
+
+ m_init(0.5*CLHEP::mm); // default road width = 0.5 CLHEP::mm
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: METHOD m_init(const double & r_road_width) ::
+//:::::::::::::::::::::::::::::::::::::::::::::::::::
+
+void StraightPatRec::m_init(const double & r_road_width) {
+
+
+//coverity
+ m_b_x1_err=0.0;
+ m_b_x2_err=0.0;
+ m_b_x2=0.0;
+ m_m_x1_err=0.0;
+ m_m_x2_err=0.0;
+ m_m_x2=0.0;
+
+
+//:::::::::::::::
+//:: VARIABLES ::
+//:::::::::::::::
+
+ Amg::Vector3D null_vec(0.0, 0.0, 0.0); // auxiliary 0 vector
+
+//::::::::::::::::::
+//:: SET TIME-OUT ::
+//::::::::::::::::::
+
+ m_time_out = 10.0;
+
+//::::::::::::::::::::::::::::
+//:: SET THE MAXIMUM RADIUS ::
+//::::::::::::::::::::::::::::
+
+ m_r_max = 15.0;
+
+//::::::::::::::::::::::::
+//:: SET THE ROAD WIDTH ::
+//::::::::::::::::::::::::
+
+ m_road_width = r_road_width;
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: INITIALIZE PRIVATE VARIABLES WHICH ARE ACCESSIBLE BY METHODS ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+ m_nb_track_hits = 0;
+ m_chi2 = 0.0;
+ m_track = MTStraightLine(null_vec, null_vec, null_vec, null_vec);
+
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+//:: SET THE TRACK IN THE x1-x3 PLANE (LATER VERSION MAY ALLOW TO SET ::
+//:: USER-DEFINED VALUES) ::
+//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+ m_m_x1 = 1.0; m_b_x1 = 0.0;
+
+ m_fix_selection = false;
+
+ return;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::
+//:: METHOD tangent ::
+//::::::::::::::::::::
+
+MTStraightLine StraightPatRec::tangent(
+ const Amg::Vector3D & r_w1,
+ const double & r_r1, const double & r_sigma12,
+ const Amg::Vector3D & r_w2, const double & r_r2,
+ const double & r_sigma22, const int & r_case) const {
+
+//:::::::::::::::
+//:: VARIABLES ::
+//:::::::::::::::
+
+ double sinalpha; // auxiliary sinus of an angle
+ double cosalpha; // auxiliary sinus of an angle
+ Amg::Vector3D e2prime, e3prime; // auxiliary direction vectors
+ MTStraightLine tang; // tangent to drift circles of a hit pair
+ Amg::Vector3D p1(0.0, 0.0, 0.0), p2(0.0, 0.0, 0.0); // hit points defining a tangent
+ Amg::Vector3D null_vec(0.0, 0.0, 0.0); // auxiliary 0 vector
+ double mx1, bx1, mx2, bx2; // auxiliary track parameters
+
+//::::::::::::::::::::::::::::::::::::::::::::
+//:: CHECK WHETHER THE SELECTED CASE EXISTS ::
+//::::::::::::::::::::::::::::::::::::::::::::
+
+ if (r_case<1 || r_case>4) {
+ cerr << endl
+ << "Class StraightPatRec, "
+ << "method tangent: ERROR!\n"
+ << "Illegal case " << r_case << ", must be 1,2,3, or 4."
+ << endl;
+ exit(1);
+ }
+
+//:::::::::::::::::::::::::::::::::::::
+//:: CALCULATE THE REQUESTED TANGENT ::
+//:::::::::::::::::::::::::::::::::::::
+
+// local coordinate axis vectors //
+ e3prime = (r_w2-r_w1).unit();
+ e2prime = Amg::Vector3D(0.0, e3prime.z(), -e3prime.y());
+
+// case 1 and 2 //
+ if (r_case==1 || r_case==2) {
+
+ sinalpha = fabs(r_r2-r_r1)/(r_w2-r_w1).mag();
+ cosalpha = sqrt(1.0-sinalpha*sinalpha);
+
+ // case 1 //
+ if (r_case==1) {
+ p1 = r_w1+((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r1*
+ (cosalpha*e2prime-sinalpha*e3prime);
+ p2 = r_w2+((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r2*
+ (cosalpha*e2prime-sinalpha*e3prime);
+ }
+
+ // case 2 //
+ if (r_case==2) {
+ p1 = r_w1-((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r1*
+ (cosalpha*e2prime+sinalpha*e3prime);
+ p2 = r_w2-((1 && r_r2>=r_r1)-(1 && r_r2<r_r1))*r_r2*
+ (cosalpha*e2prime+sinalpha*e3prime);
+ }
+
+
+ }
+
+// case 3 and 4 //
+ if (r_case==3 || r_case==4) {
+
+ sinalpha = (r_r1+r_r2)/(r_w2-r_w1).mag();
+ cosalpha = sqrt(1.0-sinalpha*sinalpha);
+
+ // case 3 //
+ if (r_case==3) {
+// p1 = r_w1+r_r1*(cosalpha*e2prime-sinalpha*e3prime);
+// p2 = r_w2-r_r2*(cosalpha*e2prime-sinalpha*e3prime);
+ p1 = r_w1+r_r1*(cosalpha*e2prime+sinalpha*e3prime);
+ p2 = r_w2-r_r2*(cosalpha*e2prime+sinalpha*e3prime);
+ }
+
+ // case 4 //
+ if (r_case==4) {
+ p1 = r_w1-r_r1*(cosalpha*e2prime-sinalpha*e3prime);
+ p2 = r_w2+r_r2*(cosalpha*e2prime-sinalpha*e3prime);
+ }
+
+ }
+
+// calculation of the tangent and estimation of its errors //
+ if ((p2-p1).z()!=0.0) {
+ tang = MTStraightLine(p1, p2-p1, null_vec, null_vec);
+ } else {
+ Amg::Vector3D direction(p2-p1);
+ direction[2] = (1.0e-99);
+ tang = MTStraightLine(p1, direction, null_vec, null_vec);
+ }
+ mx1 = tang.m_x1(); bx2 = tang.b_x1();
+ mx2 = tang.m_x2(); bx2 = tang.b_x2();
+ tang = MTStraightLine(mx1, bx1, mx2, bx2,1.0, 1.0,
+ sqrt(r_sigma12+r_sigma22)/fabs(p2.z()-p1.z()),
+ sqrt(r_sigma12+p1.z()*p1.z()*(r_sigma12+r_sigma22)/
+ pow(p2.z()-p1.z(), 2)));
+ // errors in mx1 and bx1 are arbitrary since they are
+ // not used at a later stage.
+
+
+ return tang;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD fitCandidate ::
+//:::::::::::::::::::::::::
+//
+
+
+MTStraightLine StraightPatRec::fitCandidate(MuonCalibSegment & r_segment,
+ const std::vector<unsigned int> & r_selection,
+ const MTStraightLine & cand_line) const {
+
+///////////////
+// VARIABLES //
+///////////////
+ Amg::Vector3D null(0.0, 0.0, 0.0);
+ Amg::Vector3D xhat(1.0, 0.0, 0.0);
+
+ int num_selected_hits(0);
+ for (unsigned int k=0; k<r_selection.size(); k++) {
+ if(!r_selection[k]) {
+ num_selected_hits++;
+ }
+ }
+ if(num_selected_hits<3) return cand_line;
+
+ Amg::Vector3D init_dir(r_segment.direction());
+ Amg::Vector3D init_pos(r_segment.position());
+////////////////////////////////////////////
+// SET THE CORRECT SIGNED TRACK DISTANCES //
+////////////////////////////////////////////
+
+ r_segment.set(r_segment.chi2(), cand_line.positionVector(), cand_line.directionVector());
+
+ for (unsigned int k=0; k<r_segment.mdtHitsOnTrack(); k++) {
+ MTStraightLine aux_line(r_segment.mdtHOT()[k]->localPosition(),
+ xhat, null, null);
+ r_segment.mdtHOT()[k]->setDistanceToTrack(
+ cand_line.signDistFrom(aux_line),
+ r_segment.mdtHOT()[k]->sigmaDriftRadius());
+ }
+
+//////////////////////////////////////////////////////////////////////////////////////////
+//Analitical fit of linearized measurements (transformation to reference frame is used) //
+//////////////////////////////////////////////////////////////////////////////////////////
+ unsigned int NLC=2; //number of fit parameters
+
+//Normalization on z-direction
+ Amg::Vector3D dir_norm(r_segment.direction().x()/r_segment.direction().z(),
+ r_segment.direction().y()/r_segment.direction().z(), 1.0);
+
+ Amg::Vector3D refTransl(r_segment.position());
+ refTransl[0] = 0;
+ Amg::Vector3D refDir(dir_norm);
+ refDir[0] = 0;
+ refDir = refDir.unit();
+
+// Rotation Matrix for Toewr_to_track transformation
+
+ double rotAngle = Amg::angle(Amg::Vector3D::UnitZ(),refDir)*(refDir.y())/fabs(refDir.y());
+ Amg::AngleAxis3D RotMatr(rotAngle, Amg::Vector3D::UnitX() ); // Matrix for Track reference transformation
+ Amg::AngleAxis3D RotMatr_inv(-rotAngle, Amg::Vector3D::UnitX() ); // inverse
+
+ Amg::VectorX alpha = Amg::VectorX(NLC);
+ alpha.setZero();
+ Amg::VectorX betha = Amg::VectorX(NLC);
+ betha.setZero();
+ Amg::MatrixX Gamma = Amg::MatrixX(NLC,NLC);
+ Gamma.setZero();
+// vector to store hit positions
+ vector<Amg::Vector3D > hit_position_track;
+
+ for(unsigned int l=0;l<r_segment.mdtHitsOnTrack();l++) {
+// Transformation to track reference frame
+ Amg::Vector3D WirPosLocal = r_segment.mdtHOT()[l]->localPosition();
+ Amg::Vector3D WirPosTrack = WirPosLocal-refTransl;
+ WirPosTrack = RotMatr*WirPosTrack;
+ double signedDrifRadius = fabs(r_segment.mdtHOT()[l]->driftRadius())*
+ (r_segment.mdtHOT()[l]->signedDistanceToTrack()/
+ fabs(r_segment.mdtHOT()[l]->signedDistanceToTrack()));
+ Amg::Vector3D HitPosTrack = WirPosTrack;
+ HitPosTrack[0] = HitPosTrack.y() + signedDrifRadius;
+ HitPosTrack[1] = r_segment.mdtHOT()[l]->sigma2DriftRadius(); //trick to store hit resolution
+
+ hit_position_track.push_back(HitPosTrack);
+
+ Amg::VectorX delta = Amg::VectorX(NLC);
+ delta.setZero();
+ delta[0] = 1.0;
+ delta[1] = HitPosTrack.z();
+
+ if(!r_selection[l]) {
+ double weight=1.0/(r_segment.mdtHOT()[l]->sigma2DriftRadius());
+ Gamma += weight*delta*delta.transpose();
+ betha += weight*HitPosTrack.y()*delta;
+ }
+ }
+
+// solution of linear system of equations
+ Gamma = Gamma.inverse();
+ alpha = Gamma*betha;
+
+ double refit_chi2(0.0);
+ for(unsigned int j=0;j<hit_position_track.size();j++) {
+ double res = (alpha[0] + alpha[1]*hit_position_track[j].z()
+ - hit_position_track[j].y());
+ if(!r_selection[j]) refit_chi2 += res*res/hit_position_track[j].x();
+ }
+
+// Backward transformation
+ Amg::Vector3D seg_dir(0.0,alpha[1],1.0);
+ seg_dir = RotMatr_inv*seg_dir;
+ seg_dir[1] = seg_dir.y()/seg_dir.z();
+ seg_dir[0] = init_dir.x()/init_dir.z();
+ seg_dir[2] = 1.0;
+
+ Amg::Vector3D seg_pos(0.0,alpha[0],0.0);
+ seg_pos = RotMatr_inv*seg_pos;
+ seg_pos = seg_pos + refTransl;
+
+ seg_pos[1] = seg_pos.y()+seg_dir.y()*(init_pos.z()-seg_pos.z());
+ seg_pos[0] = init_pos.x();
+ seg_pos[2] = init_pos.z();
+
+
+// Rewriting segment
+ r_segment.set(refit_chi2/static_cast<double>(num_selected_hits-2), seg_pos, seg_dir);
+
+ MTStraightLine aux_line(seg_pos, seg_dir,
+ Amg::Vector3D(0.0, 0.0, 0.0), Amg::Vector3D(0.0, 0.0, 0.0));
+
+ for (unsigned int k=0; k<r_segment.mdtHitsOnTrack(); k++) {
+ MTStraightLine aux_t(r_segment.mdtHOT()[k]->localPosition(), xhat, null, null);
+ r_segment.mdtHOT()[k]->setDistanceToTrack(aux_line.signDistFrom(aux_t),
+ r_segment.mdtHOT()[k]->sigmaDriftRadius());
+ }
+
+ return aux_line;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD roadWidth ::
+//::::::::::::::::::::::
+
+double StraightPatRec::roadWidth(void) const {
+
+ return m_road_width;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::
+//:: METHOD numberOfTrackHits ::
+//::::::::::::::::::::::::::::::
+
+unsigned int StraightPatRec::numberOfTrackHits(void) const {
+
+ return m_nb_track_hits;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::
+//:: METHOD trackHits ::
+//::::::::::::::::::::::
+
+const vector<const MdtCalibHitBase*> & StraightPatRec::trackHits(
+ void) const {
+
+ return m_track_hits;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::
+//:: METHOD chi2 ::
+//:::::::::::::::::
+
+double StraightPatRec::chi2(void) const {
+
+ if (m_nb_track_hits == 0) {
+ return -1;
+ }
+ return m_chi2;
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::::::::::::::::::::
+//:: METHOD chi2PerDegreesOfFreedom ::
+//::::::::::::::::::::::::::::::::::::
+
+double StraightPatRec::chi2PerDegreesOfFreedom(void) const{
+
+ if (m_nb_track_hits<3) {
+ return -1;
+ }
+ return chi2()/static_cast<double>(numberOfTrackHits()-2);
+
+
+}
+
+//*****************************************************************************
+
+//::::::::::::::::::
+//:: METHOD track ::
+//::::::::::::::::::
+
+MTStraightLine StraightPatRec::track(void) const {
+
+ return m_track;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::::
+//:: METHOD setRoadWidth ::
+//:::::::::::::::::::::::::
+
+void StraightPatRec::setRoadWidth(
+ const double & r_road_width) {
+
+ m_road_width = fabs(r_road_width);
+ return;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD setTimeOut ::
+//:::::::::::::::::::::::
+
+void StraightPatRec::setTimeOut(const double & time_out) {
+
+ m_time_out = time_out;
+ return;
+
+}
+
+void StraightPatRec::setFixSelection(bool fix_sel) {
+ m_fix_selection = fix_sel;
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::
+//:: METHOD fit(.) ::
+//:::::::::::::::::::
+
+bool StraightPatRec::fit(MuonCalibSegment & r_segment) const {
+
+// select all hits //
+ HitSelection selection(r_segment.mdtHitsOnTrack(), 0);
+
+// call the other fit function //
+ return fit(r_segment, selection);
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::
+//:: METHOD fit(.,.) ::
+//:::::::::::::::::::::
+bool StraightPatRec::fit(MuonCalibSegment & r_segment,
+
+ HitSelection r_selection) const {
+
+ return fitCallByReference(r_segment, r_selection);
+
+}
+
+//*****************************************************************************
+
+bool StraightPatRec::fitCallByReference(MuonCalibSegment & r_segment,
+
+ HitSelection &r_selection) const {
+
+///////////////
+// VARIABLES //
+///////////////
+
+ time_t start,end; // start and end time (needed for time-out)
+ time(&start);
+ double diff; // difference of start and end time (needed for time-out)
+ Combination combination;
+
+ vector<unsigned int> hit_index; // hit indices for given selection
+ unsigned int try_nb_hits; // try this given number of hits for the
+ // segment reconstruction
+
+ vector<const MdtCalibHitBase*> selected_hits;
+ vector<unsigned int> selected_hits_index;
+
+ Amg::Vector3D w_min, w_max; // wire with the minimum local z coordinate,
+ // wire with the maximum local z coordinate
+ double r_min, r_max; // corresponding drift CLHEP::radii
+ double sigma2_min, sigma2_max; // corresponding spatial resolution
+
+ unsigned int counter1; // auxiliary counter
+
+ unsigned int nb_candidates(0); // number of straight-track candidates
+
+ MTStraightLine tangents[4]; // the four tangents to the drift circles of a
+ // hit pair
+ MTStraightLine aux_track; // auxiliary track
+ Amg::Vector3D null(0.0, 0.0, 0.0);
+ Amg::Vector3D xhat(1.0, 0.0, 0.0);
+
+ MTStraightLine initial_track(r_segment.position(),
+ r_segment.direction(), null, null);
+ // initial track stored in the segment
+ Amg::Vector3D aux_pos, aux_dir; // auxiliary position and direction vectors
+
+///////////
+// RESET //
+///////////
+
+ m_chi2 = -1.0;
+
+ if (r_segment.mdtHitsOnTrack() != r_selection.size()) {
+ cerr << "\n"
+ << "Class StraightPatRec, "
+ << "METHOD fit(., .): WARNING!\n"
+ << "Vector with selected hits unequal to the number "
+ << "of hits on the segment!\n"
+ << "The user selection will be ignored!\n";
+ r_selection.clear();
+ r_selection.assign(r_segment.hitsOnTrack(), 0);
+ }
+
+// Filling selected_hits vector
+ for (unsigned int k=0; k<r_segment.mdtHitsOnTrack(); k++) {
+ if(!r_selection[k]) {
+ selected_hits.push_back(r_segment.mdtHOT()[k]);
+ selected_hits_index.push_back(k);
+ }
+ }
+
+///////////////////////////////////////////////////////////////////////////
+// RETURN, IF THERE ARE LESS THAN 3 HITS, I.E. THE TRACK IS NOT UNIQUELY //
+// DEFINED BY THE HITS. //
+///////////////////////////////////////////////////////////////////////////
+
+ if (selected_hits.size()<3) {
+// cerr << "\n"
+// << "Class StraightPatRec, "
+// << "METHOD fit(., .): WARNING!\n"
+// << "Too few hits for the track reconstructions!\n";
+ return false;
+ }
+
+//////////////////////////////////////////////////////////////////////////
+// FIX POTENTIAL SECOND-COORDINATE PROBLEM OF THE INITIAL TRACK SEGMENT //
+//////////////////////////////////////////////////////////////////////////
+
+ if (r_segment.direction().z()==0.0) {
+ Amg::Vector3D dir(r_segment.direction().x(),
+ r_segment.direction().y(),
+ 1.0);
+ initial_track = MTStraightLine(r_segment.position(), dir,
+ null, null);
+ }
+
+//////////////////////////
+// TRACK RECONSTRUCTION //
+//////////////////////////
+// try to find a segment with as many hits as selected initially //
+
+ HitSelection final_selection(r_selection);
+
+ try_nb_hits = selected_hits.size();
+//=============================================================================
+ while(nb_candidates==0) {
+//-----------------------------------------------------------------------------
+// Return in case of large amount of fake hits
+if(try_nb_hits<5 && (selected_hits.size()-try_nb_hits)>2) return false;
+// reset //
+ double best_chi2ndf = -1.0;
+ MTStraightLine best_aux_line(initial_track);
+
+// loop over all combinations //
+ combination.setNewParameters(selected_hits.size(), try_nb_hits);
+
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+ for (unsigned int cb=0; cb<combination.numberOfCombinations(); cb++) {
+//.............................................................................
+
+ // time-out //
+ time (&end);
+ diff = difftime (end,start);
+ if (diff>m_time_out) {
+// cerr << endl
+// << "Class StraightPatRec: "
+// << "time-out for track finding after "
+// << m_time_out
+// << " seconds!\n";
+ return false;
+ }
+
+ // get the present hit combination //
+ if (cb==0) {
+ combination.currentCombination(hit_index);
+ } else {
+ combination.nextCombination(hit_index);
+ }
+
+// Current selection of hits
+ HitSelection tmp_selection(r_selection);
+ tmp_selection.assign(tmp_selection.size(), 1);
+ for (unsigned int k=0; k<try_nb_hits; k++) {
+ tmp_selection[selected_hits_index[hit_index[k]-1]]=0;
+ }
+
+ // investigate the current combination //
+ // find the wires with minimum and maximum local z //
+ w_min = selected_hits[hit_index[0]-1]->localPosition();
+ w_max = selected_hits[hit_index[0]-1]->localPosition();
+ r_min = selected_hits[hit_index[0]-1]->driftRadius();
+ r_max = selected_hits[hit_index[0]-1]->driftRadius();
+ sigma2_min = selected_hits[hit_index[0]-1]->sigma2DriftRadius();
+ sigma2_max = selected_hits[hit_index[0]-1]->sigma2DriftRadius();
+ for (unsigned int k=1; k<try_nb_hits; k++) {
+ if (selected_hits[hit_index[k]-1]->localPosition().z()<
+ w_min.z()) {
+ w_min = selected_hits[hit_index[k]-1]->localPosition();
+ r_min = selected_hits[hit_index[k]-1]->driftRadius();
+ sigma2_min = selected_hits[hit_index[k]-1
+ ]->sigma2DriftRadius();
+ }
+ if (selected_hits[hit_index[k]-1]->localPosition().z()>
+ w_max.z()) {
+ w_max = selected_hits[hit_index[k]-1]->localPosition();
+ r_max = selected_hits[hit_index[k]-1]->driftRadius();
+ sigma2_max = selected_hits[hit_index[k]-1
+ ]->sigma2DriftRadius();
+ }
+ }
+
+ // set the spatial resolution to 0.1 CLHEP::mm if it is 0 //
+ if (sigma2_min==0) {
+ sigma2_min=0.1;
+ }
+ if (sigma2_max==0) {
+ sigma2_max=0.1;
+ }
+
+ // get the four segment candidates tangential to the outermost hits //
+ for (unsigned int r_case=1; r_case<5; r_case++) {
+ tangents[r_case-1] = tangent(w_min, r_min, sigma2_min,
+ w_max, r_max, sigma2_max, r_case);
+ }
+
+ // determine additional track points within the road width around //
+ // the four tangents and determine a segment //
+ for (unsigned int r_case=1; r_case<5; r_case++) {
+ counter1 = 0;
+ for (unsigned int n=0; n<try_nb_hits; n++) {
+ MTStraightLine aux_line(
+ selected_hits[hit_index[n]-1]->localPosition(),
+ xhat, null, null);
+ if (fabs(fabs(tangents[r_case-1].signDistFrom(
+ aux_line))-
+ selected_hits[hit_index[n]-1]->driftRadius())
+ <=m_road_width) {
+ counter1 = counter1+1;
+ }
+ }
+
+// Used to fix initial configuration of hits
+ if(m_fix_selection) counter1=try_nb_hits;
+// perform a straight line fit for the candidate //
+ if (counter1==try_nb_hits) {
+ nb_candidates++;
+ aux_track = fitCandidate(r_segment, tmp_selection,
+ tangents[r_case-1]);
+
+ if ((best_chi2ndf==-1) || (best_chi2ndf>r_segment.chi2())) {
+ best_chi2ndf = r_segment.chi2();
+ best_aux_line = aux_track;
+ final_selection.assign(final_selection.size(),1);
+ for (unsigned int j=0; j<try_nb_hits; j++) {
+ final_selection[selected_hits_index[hit_index[j]-1]]=0;
+ }
+
+ }
+ }
+ }
+//.............................................................................
+ }
+//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
+
+
+//////////////////////////////////////////////////
+// MAKE THE FINAL STRAIGHT SEGMENT, IF POSSIBLE //
+//////////////////////////////////////////////////
+
+ if (nb_candidates>0) {
+
+ // store track hits, rewrite hit selection //
+ m_track_hits.clear();
+ for (unsigned int k=0; k<r_selection.size(); k++) {
+ r_selection[k] = final_selection[k];
+ if(!final_selection[k]) {
+ m_track_hits.push_back(r_segment.mdtHOT()[k]);
+ }
+ }
+
+// Final refit //
+ m_track = fitCandidate(r_segment, final_selection,
+ best_aux_line);
+
+ m_nb_track_hits = m_track_hits.size();
+ m_chi2 = r_segment.chi2()*(static_cast<double>(m_nb_track_hits-2));
+
+ if(m_refine_segment) {
+ r_segment.refineMdtSelection(r_selection);
+ }
+
+ } else {
+ try_nb_hits = try_nb_hits-1;
+ if (try_nb_hits<3) {
+ return false;
+ }
+ }
+
+//-----------------------------------------------------------------------------
+ }
+//=============================================================================
+
+ if (std::isnan(m_chi2)) {
+ m_chi2 = -1.0;
+ return false;
+ }
+
+ return true;
+
+}
+
+//*****************************************************************************
+
+//:::::::::::::::::::::::
+//:: METHOD printLevel ::
+//:::::::::::::::::::::::
+
+void StraightPatRec::printLevel(int level) {
+
+ cerr << "\n"
+ << "Class StraightPatRec, "
+ << "method printLevel: WARNING!\n"
+ << "Print level " << level << " is ignored.\n";
+ return;
+
+}