DecorrelateMatrix.h

Go to the documentation of this file.

#include "OscLib/OscCalcSterile.h"

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Sample.h"
#include "CAFAna/Core/Registry.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "CAFAna/Cuts/QuantileCuts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Analysis/Exposures.h"

#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"

#include "NuXAna/Cuts/NusCuts18.h"
#include "NuXAna/Cuts/NusCuts20.h"
#include "NuXAna/Prediction/NDPredictionSterile.h"
#include "NuXAna/Prediction/FDPredictionSterile.h"
#include "NuXAna/Analysis/SystMaker.h"

#include "TFile.h"
#include "TKey.h"


std::vector<double> DecorrelateND(TMatrixD* inmat, int ndbins, int fdbins){

        TMatrixD matrix = *inmat;

        // Create and populate submatrices
  TMatrixD s11(ndbins, ndbins);
  TMatrixD s12(fdbins, ndbins);
  TMatrixD s21(ndbins, fdbins);
  TMatrixD s22(fdbins, fdbins);

  for (int i = 0; i < ndbins; ++i) {
    for (int j = 0; j < ndbins; ++j) {
      s11(i,j) = matrix(i,j);
    }
  }

  for (int i = 0; i < ndbins; ++i) {
    for (int j = 0; j < fdbins; ++j) {
      s12(j,i) = matrix(i,ndbins+j);
      s21(i,j) = matrix(i,ndbins+j);
    }
  }

  for (int i = 0; i < fdbins; ++i) {
    for (int j = 0; j < fdbins; ++j) {
      s22(i,j) = matrix(ndbins+i, ndbins+j);
    }
  }

  // Decorrelate
  TMatrixD tmp(ndbins, ndbins);
  s22.Invert();
  s21 *= s22;
  tmp.Mult(s21, s12);
  s11 -= tmp;

  std::vector<double> ret(ndbins);
  for (int i = 0; i < ndbins; ++i) {
    ret[i] = s11(i, i);
                //std::cout << "nd error bin " << i << ret[i] << std::endl;
  }

        return ret;
}

std::vector<double> DecorrelateFD(TMatrixD* inmat, int ndbins, int fdbins){

        TMatrixD matrix = *inmat;

  TMatrixD s11(fdbins, fdbins);
  TMatrixD s12(ndbins, fdbins);
  TMatrixD s21(fdbins, ndbins);
  TMatrixD s22(ndbins, ndbins);

  for (int i = 0; i < ndbins; ++i) {
    for (int j = 0; j < ndbins; ++j) {
      s22(i,j) = matrix(i,j);
    }
  }

  for (int i = 0; i < ndbins; ++i) {
    for (int j = 0; j < fdbins; ++j) {
      s12(i,j) = matrix(i,ndbins+j);
      s21(j,i) = matrix(i,ndbins+j);
    }
  }

  for (int i = 0; i < fdbins; ++i) {
    for (int j = 0; j < fdbins; ++j) {
      s11(i,j) = matrix(ndbins+i, ndbins+j);
    }
  }

  TMatrixD tmp(fdbins, fdbins);
  s22.Invert();
  s21 *= s22;
  tmp.Mult(s21, s12);
  s11 -= tmp;

  std::vector<double> ret(fdbins);
  for (int i = 0; i < fdbins; ++i) {
    ret[i] = s11(i, i);
                //std::cout << "fd error bin " << i << ret[i] << std::endl;
  }

        return ret;

}

std::vector<double> GetDecorrelatedUncertainty(TMatrixD* inmat, 
                                                                                                                                                                                         int ndbins, int fdbins){

        std::vector<double> decorrnd = DecorrelateND(inmat, ndbins, fdbins);
        std::vector<double> decorrfd = DecorrelateFD(inmat, ndbins, fdbins);

        std::vector<double> decorrerr;
        for (size_t i = 0; i < decorrnd.size(); ++i){
                decorrerr.push_back(std::sqrt(decorrnd.at(i)));
        }
        for (size_t i = 0; i < decorrfd.size(); ++i){
                decorrerr.push_back(std::sqrt(decorrfd.at(i)));
        }

        return decorrerr;
}