← Back to branch summary

~mwinter4/maus/ckov-update

« back to all changes in this revision

Viewing changes to src/legacy/Optics/CovarianceMatrix.cc

  • Committer: Chris Rogers
  • Date: 2012-11-23 11:19:42 UTC
  • mfrom: (659.1.50 release-candidate)
  • Revision ID: chris.rogers@stfc.ac.uk-20121123111942-5wi23zn7zx396q2d
Tags: MAUS-v0.4.1
MAUS-v0.4.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/legacy/Optics/CovarianceMatrix.cc
6
6
#include "Config/MiceModule.hh"
7
7
#include "Config/ModuleConverter.hh"
8
8
#include "Interface/Differentiator.hh"
9
 
#include "Maths/PolynomialVector.hh"
 
9
#include "Maths/PolynomialMap.hh"
10
10
#include "src/legacy/Optics/PhaseSpaceVector.hh"
11
11
 
12
12
 
798
798
  out << m_vars << "\n ************* " << std::endl;
799
799
  for(unsigned int i=0; i<m_higherMoments.size(); i++)
800
800
    out << *(m_higherMoments[i]) << "\n ************* " << std::endl;*/
801
 
  for(int i=1; i<MaxOrder()+1; i++) 
 
801
  for(int i = 0; i < MaxOrder(); ++i) 
802
802
  {
803
803
    int kvec_front = 1;
804
 
    for(unsigned int j=MAUS::PolynomialVector::NumberOfPolynomialCoefficients(6, i); j<MAUS::PolynomialVector::NumberOfPolynomialCoefficients(6, i+1); j++) 
805
 
    {
806
 
      std::vector<int> kvec = MAUS::PolynomialVector::IndexByVector(j, 6);
 
804
    for(size_t j = MAUS::PolynomialMap::NumberOfPolynomialCoefficients(6, i);
 
805
        j < MAUS::PolynomialMap::NumberOfPolynomialCoefficients(6, i+1);
 
806
        ++j) {
 
807
      std::vector<int> kvec = MAUS::PolynomialMap::IndexByVector(j, 6);
807
808
      if(kvec.front() != kvec_front) {std::cout << "\n"; kvec_front = kvec.front();}
808
809
      for(unsigned int k=0; k<kvec.size()-1; k++) std::cout << kvec[k] << ".";
809
810
      double mom = GetMoment(kvec);
822
823
  return heap.Print(out);
823
824
}
824
825
 
825
 
MAUS::PolynomialVector MomentHeap::Weighting(MomentHeap in, MomentHeap target, int order)
 
826
MAUS::PolynomialMap MomentHeap::Weighting(MomentHeap in, MomentHeap target, int order)
826
827
{
827
828
  size_t dimension = 6;
828
 
  size_t size      = MAUS::PolynomialVector::NumberOfPolynomialCoefficients(dimension, order+1);
 
829
  size_t size = MAUS::PolynomialMap::NumberOfPolynomialCoefficients(
 
830
      dimension, order);
829
831
  MAUS::Vector<double> u(size-1);
830
832
  MAUS::Matrix<double> M(size-1, size-1);
831
833
  for(size_t i=1; i<size; i++)
832
834
  {
833
 
    std::vector<int> index1 = MAUS::PolynomialVector::IndexByVector(i, dimension);
 
835
    std::vector<int> index1 = MAUS::PolynomialMap::IndexByVector(i, dimension);
834
836
    u(i)                  = target.GetMoment(index1) - in.GetMoment(index1);
835
837
    for(size_t j=1; j<size; j++)
836
838
    {
837
 
      std::vector<int> index2 = MAUS::PolynomialVector::IndexByVector(j, dimension);
 
839
      std::vector<int> index2 = MAUS::PolynomialMap::IndexByVector(j, dimension);
838
840
      std::vector<int> index3 = index2;
839
841
      index3.insert(index3.begin(), index1.begin(), index1.end());
840
842
      M(j,i)             = in.GetMoment(index3) - in.GetMoment(index1)*target.GetMoment(index2);
843
845
  MAUS::Vector<double> a  = inverse(M) * u;
844
846
  MAUS::Vector<double> a2 = MAUS::Vector<double>(a.size()+1, 1.);
845
847
  for(size_t i=0; i<a.size(); i++) a2(i+2) = a(i+1);
846
 
  return MAUS::PolynomialVector(dimension, transpose(a2) );
 
848
  return MAUS::PolynomialMap(dimension, transpose(a2) );
847
849
}