sndShowerTools.cxx File Reference

#include "sndShowerTools.h"
#include <vector>
#include <algorithm>
#include <utility>
#include "sndConfiguration.h"
#include "sndScifiPlane.h"
#include "sndUSPlane.h"
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"
#include "TVectorD.h"
#include "TDecompChol.h"
#include "Math/Vector3D.h"
#include "Math/Point3D.h"

Include dependency graph for sndShowerTools.cxx:

Go to the source code of this file.

Functions
std::pair< double, double >	PCA (const std::vector< double > &u, const std::vector< double > &z)

Function Documentation

PCA()

std::pair< double, double > PCA	(	const std::vector< double > &	u,
		const std::vector< double > &	z
	)

Definition at line 152 of file sndShowerTools.cxx.

{
  // Need at least 3 measurements
  assert(u.size() == z.size() && u.size() >= 2);
  const size_t n_samples = u.size();
  const int n_variables = 2;
 
  // Find the mean values
  double mean_u = 0.0, mean_z = 0.0;
  for (size_t i = 0; i < n_samples; ++i)
  {
    mean_u += u[i];
    mean_z += z[i];
  }
  mean_u /= n_samples;
  mean_z /= n_samples;
 
  // Find the covariance components
  double cov_uu = 0.0, cov_zz = 0.0, cov_uz = 0.0;
  for (size_t i = 0; i < n_samples; ++i)
  {
    double du = u[i] - mean_u;
    double dz = z[i] - mean_z;
    cov_uu += du * du;
    cov_zz += dz * dz;
    cov_uz += du * dz;
  }
 
  const double denom = (n_samples > 1) ? (n_samples - 1.0) : 1.0;
  cov_uu /= denom;
  cov_zz /= denom;
  cov_uz /= denom;
 
  // Fill the covariance matrix
  TMatrixDSym cov_matrix(n_variables);
  cov_matrix[0][0] = cov_uu;
  cov_matrix[0][1] = cov_uz;
  cov_matrix[1][0] = cov_uz;
  cov_matrix[1][1] = cov_zz;
 
  // Find the eigenvalues and eigenvectors. For now we only need the former.
  // Perhaps the eigenvectors could be used in the future, so left that part in comments.
  TMatrixDSymEigen eig_decomp(cov_matrix);
  TVectorD evals = eig_decomp.GetEigenValues();
  //TMatrixD evecs = eig_decomp.GetEigenVectors();
 
  // Make sure the order is descending, swap if needed
  if (evals[0] < evals[1])
  {
    std::swap(evals[0], evals[1]);
    //evecs[0] = eig_decomp.GetEigenVectors()[1];
    //evecs[1] = eig_decomp.GetEigenVectors()[0];
  }
 
  return {evals[0], evals[1]};
}