Newer
Older
TB_Chris / TbAlignment / src / .svn / text-base / TbMillepede.h.svn-base
  1. #pragma once
  2.  
  3. // Tb/TbEvent
  4. #include "Event/TbTrack.h"
  5.  
  6. // Local
  7. #include "TbAlignmentBase.h"
  8.  
  9. /** @class TbMillepede TbMillepede.h
  10. *
  11. * Implementation of the Millepede algorithm.
  12. *
  13. * @author Christoph Hombach
  14. * @date 2012-06-19
  15. */
  16.  
  17. class TbMillepede : public TbAlignmentBase {
  18. public:
  19. /// Constructor
  20. TbMillepede(const std::string& type, const std::string& name,
  21. const IInterface* parent);
  22. /// Destructor
  23. virtual ~TbMillepede();
  24.  
  25. virtual StatusCode initialize();
  26.  
  27. void align(std::vector<TbAlignmentTrack*>& alignmentTracks);
  28.  
  29. virtual void updateGeometry();
  30.  
  31. private:
  32. struct Equation {
  33. double rmeas;
  34. double weight;
  35. std::vector<int> indG;
  36. std::vector<double> derG;
  37. std::vector<int> indL;
  38. std::vector<double> derL;
  39. std::vector<int> indNL;
  40. std::vector<double> derNL;
  41. std::vector<double> slopes;
  42. };
  43. struct Constraint {
  44. /// Right-hand side (Lagrange multiplier)
  45. double rhs;
  46. /// Coefficients
  47. std::vector<double> coefficients;
  48. };
  49.  
  50. /// (Re-)initialise matrices and vectors.
  51. bool reset(const unsigned int nPlanes, const double startfact);
  52. /// Setup the constraint equations.
  53. void setConstraints(const unsigned int nPlanes);
  54. /// Define a single constraint equation
  55. void addConstraint(const std::vector<double>& dercs, const double rhs);
  56. /// Add the equations for one track and do the local fit.
  57. bool putTrack(LHCb::TbTrack* track, const unsigned int nPlanes);
  58. /// Store the parameters for one measurement.
  59. void addEquation(std::vector<Equation>& equations,
  60. const std::vector<double>& derlc,
  61. const std::vector<double>& dergb,
  62. const std::vector<double>& dernl,
  63. const std::vector<int>& dernli,
  64. const std::vector<double>& dernls, const double rmeas,
  65. const double sigma);
  66. /// Perform local parameters fit using the equations for one track.
  67. bool fitTrack(const std::vector<Equation>& equations,
  68. std::vector<double>& trackParams, const bool singlefit,
  69. const unsigned int iteration);
  70. /// Perform global parameters fit.
  71. bool fitGlobal();
  72. /// Print the results of the global parameters fit.
  73. bool printResults();
  74.  
  75. /// Matrix inversion and solution for global fit.
  76. int invertMatrix(std::vector<std::vector<double> >& v, std::vector<double>& b,
  77. const int n);
  78. /// Matrix inversion and solution for local fit.
  79. int invertMatrixLocal(std::vector<std::vector<double> >& v,
  80. std::vector<double>& b, const int n);
  81.  
  82. /// Return the limit in chi2 / ndof for n sigmas.
  83. double chi2Limit(const int n, const int nd) const;
  84. /// Multiply matrix and vector
  85. bool multiplyAX(const std::vector<std::vector<double> >& a,
  86. const std::vector<double>& x, std::vector<double>& y,
  87. const unsigned int n, const unsigned int m);
  88. /// Multiply matrices
  89. bool multiplyAVAt(const std::vector<std::vector<double> >& v,
  90. const std::vector<std::vector<double> >& a,
  91. std::vector<std::vector<double> >& w, const unsigned int n,
  92. const unsigned int m);
  93.  
  94. /// Number of global derivatives
  95. unsigned int m_nagb;
  96. /// Number of local derivatives
  97. unsigned int m_nalc;
  98.  
  99. /// Equations for each track
  100. std::vector<std::vector<Equation> > m_equations;
  101. /// Constraint equations
  102. std::vector<Constraint> m_constraints;
  103.  
  104. /// Flag for each global parameter whether it is fixed or not.
  105. std::vector<bool> m_fixed;
  106. /// Sigmas for each global parameter.
  107. std::vector<double> m_psigm;
  108.  
  109. std::vector<std::vector<double> > m_cgmat;
  110. std::vector<std::vector<double> > m_corrm;
  111. std::vector<std::vector<double> > m_clcmat;
  112.  
  113. std::vector<double> m_bgvec;
  114. std::vector<double> m_corrv;
  115. std::vector<double> m_diag;
  116.  
  117. /// Difference in misalignment parameters with respect to initial values.
  118. std::vector<double> m_dparm;
  119.  
  120. /// Mapping of internal numbering to geometry service planes.
  121. std::vector<unsigned int> m_millePlanes;
  122.  
  123. bool m_debug;
  124. /// Flag to switch on/off iterations in the global fit.
  125. bool m_iterate;
  126. /// Residual cut after the first iteration.
  127. double m_rescut;
  128. /// Residual cut in the first iteration.
  129. double m_rescut_init;
  130. /// Factor in chisquare / ndof cut.
  131. double m_cfactr;
  132. /// Reference value for chisquare / ndof cut factor.
  133. double m_cfactref;
  134. // Number of standard deviations for chisquare / ndof cut.
  135. int m_nstdev;
  136. /// Number of "full" iterations (with geometry updates).
  137. unsigned int m_nIterations;
  138. /// Sigmas for each degree of freedom
  139. std::vector<double> m_sigmas;
  140. /// Planes to be kept fixed
  141. std::vector<unsigned int> m_fixedPlanes;
  142. /// Flag to fix all degrees of freedom or only the translations.
  143. bool m_fix_all;
  144. };