EigenUtils.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Utilities.h" // For LogLikelihood(double, double)
4 
5 #include <iostream>
6 
7 namespace ana
8 {
9  //----------------------------------------------------------------------
10  Eigen::MatrixXd EigenMatrixXdFromTMatrixD(const TMatrixD* mat)
11  {
12  Eigen::MatrixXd ret(mat->GetNrows(), mat->GetNcols());
13  // TMatrixD doesn't appear to have a GetArray()
14  for(int i = 0; i < mat->GetNrows(); ++i){
15  for(int j = 0; j < mat->GetNcols(); ++j){
16  ret.coeffRef(i, j) = (*mat)(i, j);
17  }
18  }
19  return ret;
20  }
21 
22  //----------------------------------------------------------------------
23  TMatrixD TMatrixDFromEigenMatrixXd(const Eigen::MatrixXd& mat)
24  {
25  TMatrixD ret(mat.rows(), mat.cols());
26  // TMatrixD doesn't appear to have a GetArray()
27  for(int i = 0; i < mat.rows(); ++i){
28  for(int j = 0; j < mat.cols(); ++j){
29  ret(i, j) = mat.coeffRef(i, j);
30  }
31  }
32  return ret;
33  }
34 
35  //----------------------------------------------------------------------
36  double LogLikelihood(const Eigen::ArrayXd& ea, const Eigen::ArrayXd& oa, bool useOverflow)
37  {
38  assert(ea.size() == oa.size());
39 
40  double chi = 0;
41 
42  const int bufferBins = useOverflow ? 0 : -1;
43 
44  for(int i = 0; i < ea.size()+bufferBins; ++i){
45  chi += LogLikelihood(ea[i], oa[i]);
46  }
47 
48  return chi;
49  }
50 
51  //----------------------------------------------------------------------
52  double Chi2CovMx(const Eigen::ArrayXd& e, const Eigen::ArrayXd& o, const Eigen::MatrixXd& covmxinv)
53  {
54  assert(e.size() == covmxinv.rows()+2);
55 
56  Eigen::ArrayXd diff = e-o;
57  // Drop underflow and overflow bins
58  const Eigen::ArrayXd diffSub(Eigen::Map<Eigen::ArrayXd>(diff.data()+1, diff.size()-2));
59  // dot collapses things down to a single number
60  return diffSub.matrix().dot(covmxinv*diffSub.matrix());
61  }
62 
63  //----------------------------------------------------------------------
64  Eigen::MatrixXd SymmMxInverse(const Eigen::MatrixXd& mx)
65  {
66  // check if there are any null rows/columns.
67  // if there are, they make the matrix singular.
68  // we will remove them temporarily,
69  // invert the matrix, then put them back afterwards.
70  std::set<int> nullRows;
71  for (int row = 0; row < mx.rows(); row++)
72  {
73  bool rowIsNull = true;
74  for (int col = 0; col < mx.cols(); col++)
75  {
76  if (mx(row, col) != 0.)
77  {
78  rowIsNull = false;
79  break;
80  }
81  }
82 
83  if (rowIsNull)
84  nullRows.insert(row);
85  }
86 
87  std::cerr << " Notice: covariance matrix has " << nullRows.size() << " null rows.\n"
88  << "They will be removed before inverting and added back afterwards." << std::endl;
89 
90  // create a new matrix for inverting, skipping the null rows
91  Eigen::MatrixXd invMx(mx.rows() - nullRows.size(),
92  mx.cols() - nullRows.size());
93  unsigned int skippedRows = 0;
94  for (int row = 0; row < mx.rows(); row++)
95  {
96  if (nullRows.find(row) != nullRows.end())
97  {
98  skippedRows++;
99  continue;
100  }
101  unsigned int skippedCols = 0;
102  for (int col = 0; col < mx.cols(); col++)
103  {
104  // since we assumed the matrix is symmetric,
105  // we can just use the null rows list here
106  if (nullRows.find(col) != nullRows.end())
107  {
108  skippedCols++;
109  continue;
110  }
111 
112  invMx(col-skippedCols, row-skippedRows) = invMx(row-skippedRows, col-skippedCols) = mx(row, col);
113  }
114  }
115 
116  invMx = invMx.inverse();
117 
118  // put back the empty rows if there were any
119  if (!nullRows.empty())
120  {
121  skippedRows = 0;
122  Eigen::MatrixXd retMx(mx.rows(), mx.cols());
123  for (int row = 0; row < mx.rows(); row++)
124  {
125  if (nullRows.find(row) != nullRows.end())
126  {
127  skippedRows++;
128  continue;
129  }
130 
131  unsigned int skippedCols = skippedRows;
132  for (int col = 0; col < mx.cols(); col++)
133  {
134  if (nullRows.find(col) != nullRows.end())
135  {
136  skippedCols++;
137  continue;
138  }
139 
140  retMx(col, row) = retMx(row, col) = invMx(row-skippedRows, col-skippedCols);
141  }
142  }
143 
144  return retMx;
145  }
146 
147  return invMx;
148  }
149 }
double Chi2CovMx(const Eigen::ArrayXd &e, const Eigen::ArrayXd &o, const Eigen::MatrixXd &covmxinv)
Definition: EigenUtils.cxx:52
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TMatrixD TMatrixDFromEigenMatrixXd(const Eigen::MatrixXd &mat)
Definition: EigenUtils.cxx:23
OStream cerr
Definition: OStream.cxx:7
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
Int_t col[ntarg]
Definition: Style.C:29
const double j
Definition: BetheBloch.cxx:29
Eigen::MatrixXd SymmMxInverse(const Eigen::MatrixXd &mx)
Definition: EigenUtils.cxx:64
Eigen::MatrixXd EigenMatrixXdFromTMatrixD(const TMatrixD *mat)
Definition: EigenUtils.cxx:10
assert(nhit_max >=nhit_nbins)
Float_t e
Definition: plot.C:35
Eigen::MatrixXd mat