CorrelationMatrix.cc

Go to the documentation of this file.

#include <utl/CorrelationMatrix.h>
#include <utl/AugerException.h>
#include <utl/ErrorLogger.h>
#include <iostream>

using namespace std;
using namespace utl;

inline
int CorrelationMatrix::Index(const int iPar, const int jPar) const
{
  // calculates internal storage index for upper triangular matrix without diagonal elements
  if (iPar >= fNParameter || jPar >= fNParameter)
    throw OutOfBoundException("iPar or jPar too large");

  if (iPar<jPar)
    return ( (fNParameter-1)*(fNParameter-2) - 
             (fNParameter-1-iPar)*(fNParameter-2-iPar) )/2 +
           jPar - 1;
  else
    return ( (fNParameter-1)*(fNParameter-2) - 
             (fNParameter-1-jPar)*(fNParameter-2-jPar) )/2 +
           iPar - 1;
}

CorrelationMatrix::CorrelationMatrix(const int nPar) 
  : fNParameter(nPar), fElements(nPar*(nPar-1)/2)
{
  for (Elements::iterator it=fElements.begin();
       it!=fElements.end();++it)
    *it = 0.0;
}

void
CorrelationMatrix::Set(const int iPar, const int jPar, const double corr)
{
  if (iPar == jPar)
    throw DoesNotComputeException("Setting diagonal elements of correlation matrix makes no sense");

  if (corr < -1 || corr > 1)
  {
    ostringstream msg;
    msg << "Value " << corr << " is invalid (-1 <= correlation <= 1), convert to 0 (!)";
    ERROR(msg);
    
    fElements[Index(iPar,jPar)] = 0;
  } else
    fElements[Index(iPar,jPar)] = corr;
}

double
CorrelationMatrix::operator()(const int iPar, const int jPar)
  const
{
  if (iPar == jPar)
    return 1.0;

  return fElements[Index(iPar,jPar)];
}