MOISolution.cxx

Go to the documentation of this file.

#include <vector>
#include <cmath>

#include "TVectorD.h"
#include "TMatrixD.h"
#include "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"

#include "MessageService/MsgService.h"

#include "MOISolution.h"

const float C45 = sqrt(2.) / 2.; //Cosine of 45 degrees

CVSID("$Id: MOISolution.cxx,v 1.2 2007/02/04 06:16:55 rhatcher Exp $");

MOISolution::MOISolution() {
  fWF = 1.0;
  Zero();
}

MOISolution::MOISolution(const vector<double> &X, const vector<double> &Y,
                         const vector<double> &E, double WF) {
  fWF = WF;
  Fill(X, Y, E);
}

MOISolution::MOISolution(const vector<double> &X, const vector<double> &Y,
                         const vector<double> &Z, const vector<double> &E,
                         double WF) {
  fWF = WF;
  Fill(X, Y, Z, E);
}

void MOISolution::Zero() {
  EigenValues.clear();
  EigenVectors.clear();
}

void MOISolution::Fill(const vector<double> &X, const vector<double> &Y,
                       const vector<double> &E) {
  Zero();

  if (X.size() != Y.size() || X.size() != E.size() ||
      X.size() == 0) {
    MSG("MOI",Msg::kWarning) << "Mismatch size of vectors input to MOISolution" << endl;
    return;
  }
  const unsigned int NPoints = X.size();
  MSG("MOI",Msg::kVerbose) << "MOISolution::Fill 2D for " << NPoints << " points" << endl;

  double Cent[2] = {0.};

  double ESum = 0.;
  double Ethis = 0.;
  for (unsigned int i=0; i<NPoints; i++) {
    Ethis = pow(E[i],fWF);
    ESum += Ethis;

    Cent[0] += X[i] * Ethis;
    Cent[1] += Y[i] * Ethis;
  }
  if(ESum == 0.) return;//Why are you doing this to me?
  Cent[0] /= ESum;
  Cent[1] /= ESum;

  TMatrixDSym Moi(2); Moi.Zero();
  double Xthis, Ythis;
  for (unsigned int i=0; i<NPoints; i++) {
    Xthis = X[i] - Cent[0];
    Ythis = Y[i] - Cent[1];

    double Ethis = pow(E[i],fWF);

    Moi(0,0) += Ethis*(Ythis*Ythis);
    Moi(1,1) += Ethis*(Xthis*Xthis);
    Moi(0,1) += -1. * Ethis * Xthis * Ythis;
  }
  Moi(1,0)=Moi(0,1);
  Moi *= (1./ESum);

  //Do the Eigens and fill output maps
  TMatrixDSymEigen Eigen(Moi);
  TMatrixD EVec = Eigen.GetEigenVectors();
  TVectorD EVal = Eigen.GetEigenValues();

  //Dirty sort for two numbers
  int sorted[2] = {0,1};
  if (EVal(0)>=EVal(1)) {
    sorted[0] = 1; sorted[1] = 0;
  }

  for (int i=0;i<2;i++) {
    EigenValues.push_back(EVal(sorted[i]));
    vector<double> V;
    for(int j=0;j<2;j++) V.push_back(EVec(j,sorted[i]));
    EigenVectors.push_back(V);
  }
}

//Find the MOI EV
void MOISolution::Fill(const vector<double> &X, const vector<double> &Y,
                       const vector<double> &Z, const vector<double> &E) {
  Zero();
  if (X.size() != Y.size() || X.size() != Z.size() || 
      X.size() != E.size() || X.size() == 0) {
    MSG("MOI",Msg::kWarning) << "Mismatch size of vectors input to MOISolution" << endl;
    return;
  }
  const unsigned int NPoints = X.size();
  MSG("MOI",Msg::kVerbose) << "MOISolution::Fill 3D for " << NPoints << " points" << endl;

  double Cent[3] = {0.};

  double ESum = 0.;
  double Ethis = 0.;
  for (unsigned int i=0; i<NPoints; i++) {
    Ethis = pow(E[i],fWF);
    ESum += Ethis;

    Cent[0] += X[i] * Ethis;
    Cent[1] += Y[i] * Ethis;
    Cent[2] += Z[i] * Ethis;
  }
  if(ESum == 0.) return;//Why are you doing this to me?
  Cent[0] /= ESum;
  Cent[1] /= ESum;
  Cent[2] /= ESum;

  TMatrixDSym Moi(3); Moi.Zero();
  double Xthis, Ythis, Zthis;
  for (unsigned int i=0; i<NPoints; i++) {

    Xthis = X[i] - Cent[0];
    Ythis = Y[i] - Cent[1];
    Zthis = Z[i] - Cent[2];

    double Ethis = pow(E[i],fWF);

    Moi(0,0) += Ethis*(Ythis*Ythis+Zthis*Zthis);
    Moi(1,1) += Ethis*(Xthis*Xthis+Zthis*Zthis);
    Moi(2,2) += Ethis*(Xthis*Xthis+Ythis*Ythis);
    Moi(0,1) += -1. * Ethis * Xthis * Ythis;
    Moi(0,2) += -1. * Ethis * Xthis * Zthis;
    Moi(1,2) += -1. * Ethis * Ythis * Zthis;
  }
  Moi(2,0)=Moi(0,2); Moi(1,0)=Moi(0,1); Moi(2,1)=Moi(1,2);
  Moi *= (1./ESum);

  //Do the Eigens and fill the output
  TMatrixDSymEigen Eigen(Moi);
  TMatrixD EVec = Eigen.GetEigenVectors();
  TVectorD EVal = Eigen.GetEigenValues();

  //ugly sort for three numbers
  int sorted[3] = {0,1,2};
  if (EVal(0)<=EVal(1) && EVal(0)<=EVal(2)) {
    sorted[0] = 0;
    if(EVal(1)<EVal(2)) {sorted[1]=1; sorted[2]=2;}
    else {sorted[1]=2; sorted[2]=1;}
  }
  else
  if (EVal(1)<=EVal(0) && EVal(1)<=EVal(2)) {
    sorted[0] = 1;
    if(EVal(0)<EVal(2)) {sorted[1]=0; sorted[2]=2;}
    else {sorted[1]=2; sorted[2]=0;}
  }
  else
  if (EVal(2)<=EVal(0) && EVal(2)<=EVal(1)) {
    sorted[0] = 2;
    if(EVal(0)<EVal(1)) {sorted[1]=0; sorted[2]=1;}
    else {sorted[1]=1; sorted[2]=0;}
  }

  for (int i=0;i<3;i++) {
    EigenValues.push_back(EVal(sorted[i]));
    vector<double> V;
    for(int j=0;j<3;j++) V.push_back(EVec(j,sorted[i]));
    EigenVectors.push_back(V);
  }
}

---