genie_diag_utils.h

Go to the documentation of this file.

#pragma once

#include "3FlavorAna/Systs/JointAna2018Systs.h"

#include "TH1.h"
#include "TMatrixD.h"
#include "TMatrixDSymEigen.h"
#include "TMatrixDSym.h"
#include "TVectorD.h"
#include "Math/ProbFunc.h"

#include <fstream>
#include <iostream>
#include <cmath>
#include <algorithm>

double suppress = 1000000;

TH1D* ErrorBand(TH1D* hNom,std::vector<TH1D*> hUnivs,int nsigmas){
  TH1D* ret = (TH1D*)hNom->Clone(((std::string)hNom->GetName()+"_"+std::to_string(nsigmas)).c_str());

  for(int i = 1;i <= hNom->GetNbinsX();++i){
    double nom = hNom->GetBinContent(i);

    std::vector<double> events;
    for(unsigned int j = 0;j < hUnivs.size();++j)
      events.push_back(hUnivs[j]->GetBinContent(i));
    sort(events.begin(),events.end());

    double pos=0;
    for(unsigned int j = 0;j < events.size()-1;++j)
      if(nom >= events[j] && nom <= events[j+1]){
        pos = j+0.5;
        break;
      }

    double count_fraction = ROOT::Math::normal_cdf(nsigmas)-ROOT::Math::normal_cdf(0);

    int index = 0;
    int boundIdx = pos + count_fraction*(double)events.size();
    if(nsigmas >= 0) index = min(boundIdx, (int)events.size()-1);
    else index = max(boundIdx, 0);

    ret->SetBinContent(i,events[index]);
  }

  return ret;
}

TH1D* CombineHistograms(std::vector<TH1D*> hists){
  static int uniq = 0;
  std::string hist_name = "hist_combined"+std::to_string(uniq);

  int nbins = 0;
  for(unsigned int i = 0;i < hists.size();++i) nbins += hists[i]->GetNbinsX();

  TH1D* ret = new TH1D(hist_name.c_str(),"",nbins,0,nbins);

  int bin = 1;
  for(auto hist: hists){
    for(int i = 1;i <= hist->GetNbinsX();++i){
      ret->SetBinContent(bin,hist->GetBinContent(i));
      ++bin;
    }
  }

  ++uniq;
  return ret;
}

std::vector<TH1D*> SplitHistograms(TH1D* hist, int nhists)
{
  // get back the flavor spectra
  static int uniq2 = 0;
  int ret_bins = hist->GetNbinsX()/nhists;
  std::vector<TH1D*> hists;

  for(int i = 0; i < nhists; i++){
    std::string hist_name = "hist_split"+std::to_string(uniq2);
    TH1D* ret = new TH1D(hist_name.c_str(), "", ret_bins, 0, ret_bins);
    
    for(int bin = 1; bin <= ret_bins; bin++)
     ret->SetBinContent(bin, hist->GetBinContent(ret_bins*i + bin));

    hists.push_back(ret);
    uniq2++;
  } 

  return hists;
}

std::unique_ptr<TMatrixDSym> GetCovMx(const std::vector<TH1D*> &hists)
{
  // a bit different from the one in Utilities because I use TMatrixDSym
  if (hists.size() < 2)
    return std::unique_ptr<TMatrixDSym>(nullptr);

  int nBins = hists[0]->GetNbinsX();  // 1 and nBins are inclusive

  std::vector<double> binMeans(nBins);
  std::for_each(hists.begin(), hists.end(), 
                [&](TH1D* hist){
                  for(int bin = 1; bin <= nBins; bin++) 
                    binMeans[bin-1] += hist->GetBinContent(bin)/hists.size();
                });                             

  auto covMx = std::make_unique<TMatrixDSym>(nBins);
  std::for_each(hists.begin(), hists.end(), 
                [&](TH1D* hist){
                  for(int bin1 = 1; bin1 <= nBins; bin1++){
                    double x1 = hist->GetBinContent(bin1);
                    for(int bin2 = bin1; bin2 <= nBins; bin2++){
                      double x2 = hist->GetBinContent(bin2);
                      (*covMx)(bin1-1, bin2-1) += (x1-binMeans[bin1-1])*(x2-binMeans[bin2-1])/(hists.size()-1);
                      if(bin1 != bin2) 
                        (*covMx)(bin2-1,bin1-1) = (*covMx)(bin1-1,bin2-1);
                    }
                  }
                });                             


  
  // Check if its symmetric 
  if(covMx->IsSymmetric()) std::cout << "Covariance Matrix is symmetric!" << std::endl;
  assert(covMx->IsSymmetric());
 
  return covMx;
}

void CrossCheckDiag(TMatrixDSym covMx, TMatrixD V, 
                    TMatrixD eigenvectors, TVectorD eigenvalues, double tolerance = 1e-14)
{
  // Number of cross-checks after diagonalizing
  TMatrixD C = eigenvectors*V;
  TMatrixD I(TMatrixD::kUnit, V);
  std::cout << "Checking for Orthogonality.." << std::endl;
  VerifyMatrixIdentity(C, I, 1, tolerance);             // orthogonality

  std::cout << "Verifying Eigenvectors.." << std::endl;
  for(int i = 0; i < eigenvalues.GetNrows(); i++){
    TVectorD left = eigenvectors[i];
    left *= covMx;
    TVectorD right = eigenvectors[i];
    right *= eigenvalues[i];
    std::cout << "Component " << i << std::endl;
    VerifyVectorIdentity(left, right, 1, tolerance);    // eigenvector relations
  }

  std::cout << "Checking if Covariance Matrix Trace = Sum of Eigenvalues.." << std::endl;
  double trace = 0;
  double eigensum = 0;
  for(int i = 0; i < (covMx).GetNrows(); i++) trace += (covMx)(i,i);
  for(int i = 0; i < eigenvalues.GetNrows(); i++) eigensum += eigenvalues[i];
  std::cout << "covariance trace : " << sqrt(trace) << " eigenvalue sum : " << sqrt(eigensum) << std::endl;
  std::cout << std::abs(trace-eigensum) << std::endl;
  //assert(std::abs(trace-eigensum) < tolerance);         // coverage of variances
  std::cout << "Continuing.." << std::endl;

}