Belle II Software  release-05-02-19
PrecisionMatrixUtil.h
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Oliver Frost *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #pragma once
11 
12 #include <tracking/trackFindingCDC/numerics/JacobianMatrixUtil.h>
13 
14 #include <tracking/trackFindingCDC/numerics/PrecisionMatrix.h>
15 #include <tracking/trackFindingCDC/numerics/JacobianMatrix.h>
16 #include <tracking/trackFindingCDC/numerics/ParameterVector.h>
17 
18 #include <tracking/trackFindingCDC/numerics/EigenView.h>
19 
20 #include <Eigen/Core>
21 
22 #include <type_traits>
23 
24 namespace Belle2 {
29  namespace TrackFindingCDC {
30 
32  struct PrecisionMatrixUtil {
33 
35  template <int N>
36  static PrecisionMatrix<N> identity()
37  {
39  }
40 
42  template <int N>
43  static void transport(const JacobianMatrix<N, N>& ambiguity, PrecisionMatrix<N>& precision)
44  {
45  mapToEigen(precision) =
46  mapToEigen(ambiguity).transpose() * mapToEigen(precision) * mapToEigen(ambiguity);
47  }
48 
51  template <int M, int N>
52  static PrecisionMatrix<M> transported(const JacobianMatrix<N, M>& ambiguity,
53  const PrecisionMatrix<N>& precision)
54  {
55  return mapToEigen(ambiguity).transpose() * mapToEigen(precision) * mapToEigen(ambiguity);
56  }
57 
59  template <int N>
60  static void scale(const ParameterVector<N>& scales, PrecisionMatrix<N>& precision)
61  {
62  return transport(JacobianMatrixUtil::scale(1 / scales), precision);
63  }
64 
66  template <int N>
67  static PrecisionMatrix<N> scale(const ParameterVector<N>& scales,
68  const PrecisionMatrix<N>& precision)
69  {
70  return transported(JacobianMatrixUtil::scale(1 / scales), precision);
71  }
72 
75  template <int N1, int N2>
77  const PrecisionMatrix<N2>& block2)
78  {
80  mapToEigen(result) << mapToEigen(block1), Eigen::Matrix<double, N1, N2>::Zero(),
81  Eigen::Matrix<double, N2, N1>::Zero(), mapToEigen(block2);
82  return result;
83  }
84 
86  template <class APrecisionMatrix, int I = 0, int N = 0>
87  static APrecisionMatrix getSub(const PrecisionMatrix<N>& precision)
88  {
89  constexpr const int M =
90  std::remove_reference_t<decltype(mapToEigen(APrecisionMatrix()))>::RowsAtCompileTime;
91  return precision.template block<M, M>(I, I);
92  }
93 
105  template <int N>
106  static double average(const ParameterVector<N>& parameter1,
107  const PrecisionMatrix<N>& precision1,
108  const ParameterVector<N>& parameter2,
109  const PrecisionMatrix<N>& precision2,
110  ParameterVector<N>& parameter,
111  PrecisionMatrix<N>& precision)
112  {
113  const auto& ePrecision1 = mapToEigen(precision1);
114  const auto& ePrecision2 = mapToEigen(precision2);
115  auto&& ePrecision = mapToEigen(precision);
116 
117  const auto& eParameter1 = mapToEigen(parameter1);
118  const auto& eParameter2 = mapToEigen(parameter2);
119  auto&& eParameter = mapToEigen(parameter);
120 
121  ePrecision = ePrecision1 + ePrecision2;
122  eParameter = ePrecision.colPivHouseholderQr().solve(ePrecision1 * eParameter1 +
123  ePrecision2 * eParameter2);
124 
125  auto eResidual1 = eParameter1 - eParameter;
126  auto eResidual2 = eParameter2 - eParameter;
127 
128  Eigen::Matrix<double, 1, 1> chi2 = (eResidual1.transpose() * ePrecision1 * eResidual1 +
129  eResidual2.transpose() * ePrecision2 * eResidual2);
130  return chi2[0];
131  }
132 
147  template <int M, int N1, int N2>
148  static double average(const ParameterVector<N1>& parameter1,
149  const PrecisionMatrix<N1>& precision1,
150  const JacobianMatrix<N1, M>& ambiguity1,
151  const ParameterVector<N2>& parameter2,
152  const PrecisionMatrix<N2>& precision2,
153  const JacobianMatrix<N2, M>& ambiguity2,
154  ParameterVector<M>& parameter,
155  PrecisionMatrix<M>& precision)
156  {
157  const auto& eParameter1 = mapToEigen(parameter1);
158  const auto& ePrecision1 = mapToEigen(precision1);
159  const auto& eAmbiguity1 = mapToEigen(ambiguity1);
160 
161  const auto& eParameter2 = mapToEigen(parameter2);
162  const auto& ePrecision2 = mapToEigen(precision2);
163  const auto& eAmbiguity2 = mapToEigen(ambiguity2);
164 
165  auto&& ePrecision = mapToEigen(precision);
166  auto&& eParameter = mapToEigen(parameter);
167 
168  ePrecision = (eAmbiguity1.transpose() * ePrecision1 * eAmbiguity1 +
169  eAmbiguity2.transpose() * ePrecision2 * eAmbiguity2);
170 
171  eParameter = ePrecision.colPivHouseholderQr().solve(
172  eAmbiguity1.transpose() * ePrecision1 * eParameter1 +
173  eAmbiguity2.transpose() * ePrecision2 * eParameter2);
174 
175  auto eResidual1 = eParameter1 - eAmbiguity1 * eParameter;
176  auto eResidual2 = eParameter2 - eAmbiguity2 * eParameter;
177 
178  Eigen::Matrix<double, 1, 1> chi2 = (eResidual1.transpose() * ePrecision1 * eResidual1 +
179  eResidual2.transpose() * ePrecision2 * eResidual2);
180  return chi2[0];
181  }
182  };
183  }
185 }
Belle2::TrackFindingCDC::JacobianMatrixUtil::scale
static JacobianMatrix< N > scale(const ParameterVector< N > &scales)
Calculates the jacobian matrix for a scaling in each parameter.
Definition: JacobianMatrixUtil.h:76
Belle2::TrackFindingCDC::PrecisionMatrixUtil::getSub
static APrecisionMatrix getSub(const PrecisionMatrix< N > &precision)
Gets a subprecision from a precision matrix.
Definition: PrecisionMatrixUtil.h:95
Belle2::TrackFindingCDC::PrecisionMatrixUtil::scale
static void scale(const ParameterVector< N > &scales, PrecisionMatrix< N > &precision)
Scale the precision inplace by the given factors in each parameter.
Definition: PrecisionMatrixUtil.h:68
Belle2::TrackFindingCDC::PrecisionMatrixUtil::stackBlocks
static PrecisionMatrix< N1+N2 > stackBlocks(const PrecisionMatrix< N1 > &block1, const PrecisionMatrix< N2 > &block2)
Combines two precision matrices by putting them in two blocks on the diagonal of a larger matrix.
Definition: PrecisionMatrixUtil.h:84
Belle2::TrackFindingCDC::PrecisionMatrixUtil::transport
static void transport(const JacobianMatrix< N, N > &ambiguity, PrecisionMatrix< N > &precision)
Transport the precision matrix inplace with the given jacobian matrix.
Definition: PrecisionMatrixUtil.h:51
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::PrecisionMatrixUtil::identity
static PrecisionMatrix< N > identity()
Constructs an identity matrix.
Definition: PrecisionMatrixUtil.h:44
Belle2::TrackFindingCDC::PlainMatrix
A matrix implementation to be used as an interface typ through out the track finder.
Definition: PlainMatrix.h:50
Belle2::TrackFindingCDC::PrecisionMatrixUtil::average
static double average(const ParameterVector< N > &parameter1, const PrecisionMatrix< N > &precision1, const ParameterVector< N > &parameter2, const PrecisionMatrix< N > &precision2, ParameterVector< N > &parameter, PrecisionMatrix< N > &precision)
Averages two parameter vectors taking into account their respective precision.
Definition: PrecisionMatrixUtil.h:114
Belle2::TrackFindingCDC::PlainMatrix::Identity
static PlainMatrix< T, M, N > Identity()
Construct an identity matrix.
Definition: PlainMatrix.h:84
Belle2::TrackFindingCDC::PrecisionMatrixUtil::transported
static PrecisionMatrix< M > transported(const JacobianMatrix< N, M > &ambiguity, const PrecisionMatrix< N > &precision)
Return a copy of the precision matrix transported with the given back projection jacobian matrix.
Definition: PrecisionMatrixUtil.h:60