CalibrationDataGlobalEigenVariations Class Reference

#include <CalibrationDataEigenVariations.h>

Inheritance diagram for CalibrationDataGlobalEigenVariations:

Collaboration diagram for CalibrationDataGlobalEigenVariations:

Public Types
typedef std::set< size_t >	IndexSet
typedef std::set< IndexSet >	IndexSuperSet

Public Member Functions
	CalibrationDataGlobalEigenVariations (const std::string &cdipath, const std::string &tagger, const std::string &wp, const std::string &jetcollection, const std::vector< std::string > &flavours, CalibrationDataHistogramContainer *cnt, bool excludeRecommendedUncertaintySet=false)
	~CalibrationDataGlobalEigenVariations ()
TMatrixDSym	getEigenCovarianceMatrix ()
	also provide (some) access to the underlying information: covariance matrix corresponding to eigenvector variations
void	initialize (double min_variance=1.0E-6)
	carry out the eigenvector computations.
TMatrixD	getJacobianReductionMatrix (TMatrixDSym &cov)
void	excludeNamedUncertainty (const std::string &name, const std::string &flavour)
std::vector< std::string >	listNamedVariations (const std::string &flavour) const
unsigned int	getNamedVariationIndex (const std::string &name, const std::string &flavour) const
bool	getNamedVariation (const std::string &flavour, const std::string &name, TH1 *&up, TH1 *&down)
bool	getNamedVariation (unsigned int nameIndex, const std::string &flavour, TH1 *&up, TH1 *&down)
bool	getEigenvectorVariation (const std::string &flavour, unsigned int variation, TH1 *&up, TH1 *&down)
unsigned int	getNumberOfEigenVariations (const std::string &flavour)
bool	isExtrapolationVariation (unsigned int nameIndex, const std::string &flavour) const
void	mergeVariationsFrom (const size_t &index, std::string &flav)
void	mergeVariations (const IndexSet &set, std::string &flav)
void	mergeVariations (const IndexSuperSet &set, std::string &flav)
void	removeVariations (const IndexSet &set, std::string &flav)
void	removeVariations (const IndexSuperSet &set, std::string &flav)
void	excludeNamedUncertainty (const std::string &name, CalibrationDataContainer *cnt)
	exclude the source of uncertainty indicated by name from eigenvector calculations
void	removeVariations (const IndexSet &set)
	remove all variations in the given set
void	removeVariations (const IndexSuperSet &set)
	remove all variations in any of the given sets
void	mergeVariationsFrom (const size_t &index)
	merge all variations starting from the given index
void	mergeVariations (const IndexSet &set)
	merge all variations in the given set
void	mergeVariations (const IndexSuperSet &set)
	merge all variations in any of the given sets
unsigned int	getNumberOfNamedVariations () const
	retrieve the number of named variations
std::vector< std::string >	listNamedVariations () const
	list the named variations
unsigned int	getNamedVariationIndex (const std::string &name) const
	retrieve the integer index corresponding to the named variation.
unsigned int	getNumberOfEigenVariations ()
	retrieve the number of eigenvector variations
bool	getEigenvectorVariation (unsigned int variation, TH1 *&up, TH1 *&down)
	obtain the "up" and "down" variations for the given eigenvector number.
bool	getNamedVariation (const std::string &name, TH1 *&up, TH1 *&down)
	obtain the "up" and "down" variations for the named uncertainty.
bool	getNamedVariation (unsigned int nameIndex, TH1 *&up, TH1 *&down)
	obtain the "up" and "down" variations for the source uncertainty pointed to by the given index (which is assumed to correspond to the value retrieved using getNamedVariationIndex()).
bool	isExtrapolationVariation (unsigned int nameIndex) const
	flag whether the given index corresponds to an extrapolation variation
TMatrixDSym	getEigenCovarianceMatrixFromVariations ()
	covariance matrix corresponding to eigenvector variations constructed from the eigen-variation
TMatrixD	getJacobianReductionMatrix ()
	matrix to remove unecessary rows and columns from covariance
bool	EigenVectorRecomposition (const std::string &label, std::map< std::string, std::map< std::string, float > > &coefficientMap)
	Eigenvector recomposition method.
void	setVerbose (bool)

Protected Attributes
bool	m_initialized
	flag whether the initialization has been carried out
bool	m_validate
std::map< std::string, unsigned int >	m_namedIndices
	named variations
std::vector< std::pair< TH1 *, TH1 * > >	m_named
int	m_namedExtrapolation
	named variation index for the special case of extrapolation uncertainties
std::vector< std::pair< TH1 *, TH1 * > >	m_eigen
	eigenvector variations
bool	m_statVariations
	indicate whether statistical uncertainties are stored as variations
std::string	m_cdipath
	@ data members needed for eigenvector method
std::string	m_taggername
std::string	m_wp
std::string	m_jetauthor
double	m_totalvariance
double	m_capturedvariance
bool	m_verbose

Private Attributes
int	m_blockmatrixsize
std::map< std::string, Analysis::CalibrationDataHistogramContainer * >	m_histcontainers
std::set< std::string >	m_all_shared_systematics
std::set< std::string >	m_only_shared_systematics
std::map< std::string, std::vector< int > >	m_flavour_combinations
std::map< std::string, std::map< std::string, unsigned int > >	m_flav_namedIndices
std::map< std::string, std::vector< std::pair< TH1 *, TH1 * > > >	m_flav_named
std::map< std::string, int >	m_flav_namedExtrapolation
std::map< std::string, std::vector< std::pair< TH1 *, TH1 * > > >	m_flav_eigen
std::vector< std::string >	m_flavours
CalibrationDataHistogramContainer *	m_cnt
	container object containing the basic information

Detailed Description

Definition at line 145 of file CalibrationDataEigenVariations.h.

Member Typedef Documentation

IndexSet

typedef std::set<size_t> Analysis::CalibrationDataEigenVariations::IndexSet

inherited

Definition at line 29 of file CalibrationDataEigenVariations.h.

IndexSuperSet

typedef std::set<IndexSet> Analysis::CalibrationDataEigenVariations::IndexSuperSet

inherited

Definition at line 30 of file CalibrationDataEigenVariations.h.

Constructor & Destructor Documentation

CalibrationDataGlobalEigenVariations()

CalibrationDataGlobalEigenVariations::CalibrationDataGlobalEigenVariations	(	const std::string &	cdipath,
		const std::string &	tagger,
		const std::string &	wp,
		const std::string &	jetcollection,
		const std::vector< std::string > &	flavours,
		CalibrationDataHistogramContainer *	cnt,
		bool	excludeRecommendedUncertaintySet = false )

Definition at line 1290 of file CalibrationDataEigenVariations.cxx.

                                                                                                                                                                                                                                                                                                      : 
CalibrationDataEigenVariations(cdipath, tagger, wp, jetcollection, cnt, excludeRecommendedUncertaintySet, false), m_flav_eigen(), m_flavours(flavours)
//m_blockmatrixsize(1200)
{
  m_blockmatrixsize = 0; // <------- This needs to be computed based off of the dimensions of the result vectors (currently 4x1x5 for continuous WP, 1x300x1 for fixed cut WP, circa 2022)
  m_initialized = false; 
  m_namedExtrapolation = -1; 
  m_statVariations = false;
  // We want to retrieve all the containers that belong in the tagger/jetcollection/wp group
  // These are then stored, with all uncertainties extracted and stored as well, for later use in EV decomposition
 
  // For now we will open up the CDI file from within the constructor, to get at the data we need
  TString fileName(cdipath);
  TFile* f = TFile::Open(fileName.Data(), "READ");
  f->cd();
 
  // This loop extracts all the "flavour containers", i.e. all CalibrationDataHistogramContainers pertaining to the same path, but with different flavours
  // It also puts all the uncertainties for all containers into a set, which we can then later loop over, to construct the total covariance matrix
  for (const auto& flavour : m_flavours){
    std::string dir = m_taggername + "/" + m_jetauthor + "/" + m_wp  + "/" + flavour + "/" + "default_SF" ;
    TString contName(dir);
    Analysis::CalibrationDataHistogramContainer* c;
    f->GetObject(contName.Data(), c);
    if (c) {
      if (m_verbose) std::cout << "Found " << contName.Data() << std::endl;
      m_histcontainers.insert({flavour, c}); // Build the mapping between flavour and the corresponding "flavour container", i.e. the CalibrationDataHistogramContainer
      std::vector<std::string> uncs = c->listUncertainties();
      TH1* result = dynamic_cast<TH1*>(c->GetValue("result")); // let's get the size of this for later
      if (not result){
        if (m_verbose) std::cout << "Dynamic cast failed at "<<__LINE__<<"\n";
        continue;
      }
      m_blockmatrixsize+=result->GetNbinsX()*result->GetNbinsY()*result->GetNbinsZ(); // should be ~300 for fixed cut, something else for continuous
      if (m_verbose) std::cout << "m_blockmatrixsize is now " << m_blockmatrixsize << std::endl;
      for (const std::string& unc : uncs){
        if (unc.find("stat_np") != string::npos) m_statVariations = true;
        if ((unc=="result")||(unc=="comment")||(unc=="ReducedSets")||(unc=="systematics")||(unc=="statistics")||(unc=="extrapolation")||(unc=="MChadronisation")||(unc=="combined")||(unc=="extrapolation from charm")) { 
          continue; 
        } else { 
          m_all_shared_systematics.insert(unc); // Construct the set of uncertainties to get a full tally of everything in the container group
        }
      } 
      // If specified, add items recommended for exclusion from EV decomposition by the calibration group to the 'named uncertainties' list
      // This used to work on only one container, but now we do it on all four containers
      if (excludeRecommendedUncertaintySet) {
        std::vector<std::string> to_exclude = split(c->getExcludedUncertainties());
        if (to_exclude.size() == 0) {
          std::cerr << "CalibrationDataEigenVariations warning: exclusion of pre-set uncertainties list requested but no (or empty) list found" << std::endl;
        }
        for (const auto& name : to_exclude) {
          if (name == "") continue;
          excludeNamedUncertainty(name, flavour);
        }
      }
    }
  }
  
  if (m_verbose) std::cout << "\n number of shared uncertainties is " << m_all_shared_systematics.size() << std::endl;
  
  std::set<std::string>::iterator it = m_all_shared_systematics.begin();
  if (it != m_all_shared_systematics.end()){
    if (m_verbose) std::cout << "Printing out all shared uncertainties for " << tagger << "/" << jetcollection << "/" << wp << std::endl;
    if (m_verbose) std::cout << "| " << std::endl;
    while (it != m_all_shared_systematics.end()){
      if (m_verbose) std::cout << "|-- " << (*it) << std::endl;
      ++it;
    }
  } else {
    if (m_verbose) std::cout << "| no shared systematics between ";
    for (const auto& f : m_flavours){
      if (m_verbose) std::cout << f << ", ";
    } if (m_verbose) std::cout << std::endl;
  }
  // End of constructor
}

~CalibrationDataGlobalEigenVariations()

CalibrationDataGlobalEigenVariations::~CalibrationDataGlobalEigenVariations

(

)

Definition at line 1367 of file CalibrationDataEigenVariations.cxx.

{
  // delete all variation histograms owned by us
  for (const auto& flavour : m_flavours){
    for (vector<pair<TH1*, TH1*> >::iterator it = m_flav_eigen[flavour].begin(); it != m_flav_eigen[flavour].end(); ++it) {
      delete it->first;
      delete it->second;
    }
 
    for (vector<pair<TH1*, TH1*> >::iterator it = m_flav_named[flavour].begin(); it != m_flav_named[flavour].end(); ++it) {
      delete it->first;
      delete it->second;
    }
  }
}

Member Function Documentation

EigenVectorRecomposition()

bool CalibrationDataEigenVariations::EigenVectorRecomposition ( const std::string & label, std::map< std::string, std::map< std::string, float > > & coefficientMap )

inherited

Eigenvector recomposition method.

Definition at line 1113 of file CalibrationDataEigenVariations.cxx.

{
  // Calculating eigen vector recomposition coefficient map and pass to
  // user by reference. Return true if method success. Return false and
  // will not modify coefficientMap if function failed.
  //
  //     label:          flavour label
  //     coefficientMap: (reference to) coefficentMap which will be used as return value.
 
  if (! m_initialized) initialize();
 
  std::vector<TH1*> originSF_hvec;
  std::vector<std::unique_ptr<TH1>> eigenSF_hvec;
 
  // Retrieving information for calculation
  std::vector<string>fullUncList = m_cnt->listUncertainties();
  std::vector<string> uncList;
  for (unsigned int t = 0; t < fullUncList.size(); ++t) {
    // entries that should never be included
    if (fullUncList[t] == "comment" || fullUncList[t] == "result" || fullUncList[t] == "combined" || 
        fullUncList[t] == "statistics" || fullUncList[t] == "systematics" || fullUncList[t] == "MCreference" || 
        fullUncList[t] == "MChadronisation" || fullUncList[t] == "extrapolation" || fullUncList[t] == "ReducedSets" || 
        fullUncList[t] == "excluded_set") continue;
 
    // entries that can be excluded if desired
    if (m_namedIndices.find(fullUncList[t]) != m_namedIndices.end()) continue;
    
    TH1* hunc = dynamic_cast<TH1*>(m_cnt->GetValue(fullUncList[t].c_str()));
    if (not hunc){
      std::cerr<<"CalibrationDataEigenVariations::EigenVectorRecomposition: dynamic cast failed\n";
      continue;
    }
 
    Int_t nx = hunc->GetNbinsX();
    Int_t ny = hunc->GetNbinsY();
    Int_t nz = hunc->GetNbinsZ();
    Int_t bin = 0;
    bool retain = false; // Retain the histogram?
 
    // discard empty histograms
    // Read all bins without underflow&overflow
    for(Int_t binx = 1; binx <= nx; binx++)
      for(Int_t biny = 1; biny <= ny; biny++)
    for(Int_t binz = 1; binz <= nz; binz++){
      bin = hunc->GetBin(binx, biny, binz);  
      if (fabs(hunc->GetBinContent(bin)) > 1E-20){
        retain = true;
        break;
      }
    }// end hist bin for-loop
    if (!retain){
      if (m_verbose) std::cout<<"Eigenvector Recomposition: Empty uncertainty "<<fullUncList.at(t)<<" is discarded."<<std::endl;
      continue; // discard the vector
    }
 
    uncList.push_back(fullUncList.at(t));
    originSF_hvec.push_back(hunc);
  }
 
  TH1* nom = dynamic_cast<TH1*>(m_cnt->GetValue("result")); // Nominal SF hist
  if (not nom){
     if (m_verbose) std::cout<<"Eigenvector Recomposition: dynamic cast failed\n";
     return false;
  }
  int dim = nom->GetDimension();
  Int_t nx = nom->GetNbinsX();
  Int_t ny = nom->GetNbinsY();
  Int_t nz = nom->GetNbinsZ();
  Int_t nbins = nx;
  if(dim > 1) nbins *= ny;
  if(dim > 2) nbins *= nz;
  TMatrixD matSF(uncList.size(), nbins);
  Int_t col = 0; // mark the column number
  // Fill the Delta SF Matrix
  for(unsigned int i = 0; i < uncList.size(); i++){
    col = 0;
    // Loop all bins except underflow&overflow bin
    for(int binz = 1; binz <= nz; binz++)
      for(int biny = 1; biny <= ny; biny++)
    for(int binx = 1; binx <= nx; binx++){
      Int_t bin = originSF_hvec.at(i)->GetBin(binx, biny, binz);
      TMatrixDRow(matSF,i)[col] = originSF_hvec[i]->GetBinContent(bin);
      col++;
    }
  }
 
  // get eigen vectors of scale factors. Note that this is not the original eigen-vector.
  int nEigen = getNumberOfEigenVariations();
  TH1* up = nullptr;
  TH1* down = nullptr;
  for (int i = 0; i < nEigen; i++){
    if (!getEigenvectorVariation(i, up, down)){
       std::cerr<<"EigenVectorRecomposition: Error on retrieving eigenvector "<<i<<std::endl;
      return false;
    }
    // Clone first: up is owned by m_eigen, so subtracting nom in place would corrupt the cache.
    std::unique_ptr<TH1> pure(static_cast<TH1*>(up->Clone()));
    pure->SetDirectory(nullptr);
    pure->Add(nom, -1);
    eigenSF_hvec.push_back(std::move(pure));
  }
  TMatrixD matEigen(nEigen, nbins);
 
  // Fill the Eigen Matrix
  for(int i = 0; i < nEigen; i++){
    col = 0;
    // Read 300 bins without underflow&overflow
    for(int binz = 1; binz <= nz; binz++)
      for(int biny = 1; biny <= ny; biny++)
    for(int binx = 1; binx <= nx; binx++){
      Int_t bin = eigenSF_hvec.at(i)->GetBin(binx, biny, binz);
      TMatrixDRow(matEigen,i)[col] = eigenSF_hvec[i]->GetBinContent(bin);
      col++;
    }
  }
 
  eigenSF_hvec.clear();
 
  // Sanity check:
  TMatrixD matEigenOriginal = matEigen;
  TMatrixD matEigenTranspose = matEigen;
  matEigenTranspose = matEigenTranspose.T();
  TMatrixD matOriginalTimesTranspose = matEigenOriginal*matEigenTranspose;
  TMatrixD matEigenInvert = matEigenTranspose*matOriginalTimesTranspose.Invert();
  //(matEigenOriginal*matEigenInvert).DrawClone("colz"); // This should give us an identity matrix
 
  TMatrixD matCoeff = matSF*matEigenInvert;
  int nRows = matCoeff.GetNrows();
  int nCols = matCoeff.GetNcols();
  std::map<std::string, float> temp_map;
  for (int col = 0; col < nCols; col++){
    temp_map.clear();
    for(int row = 0; row < nRows; row++){
      temp_map[uncList[row]] = TMatrixDRow(matCoeff, row)[col];
    }
    coefficientMap["Eigen_"+label+"_"+std::to_string(col)] = temp_map;
  }
  
  return true;
}

excludeNamedUncertainty() [1/2]

void CalibrationDataEigenVariations::excludeNamedUncertainty ( const std::string & name, CalibrationDataContainer * cnt )

inherited

exclude the source of uncertainty indicated by name from eigenvector calculations

Definition at line 371 of file CalibrationDataEigenVariations.cxx.

{
  // Exclude the source of uncertainty identified by the given name from being used
  // in the construction of the covariance matrix to be diagonalised.
  // Notes:
  // - Some names returned by CalibrationDataContainer::listUncertainties() are not
  //   meaningful in this context, and specifying them is not allowed.
  // - Once the eigenvector diagonalisation has been carried out, this method may
  //   not be used anymore.
 
  if (m_initialized) {
    std::cerr << "CalibrationDataEigenVariations::excludeNamedUncertainty error:" << " initialization already done" << std::endl;
  } else if (name == "comment"    || name == "result"   || name == "systematics" || 
             name == "statistics" || name == "combined" || name == "extrapolation" || 
             name == "MCreference" || name == "MChadronisation" || name == "ReducedSets" || 
             name == "excluded_set" || name == "" || name == " ") // <--- these last two handle some cases that may turn up
  {
    std::cerr << "CalibrationDataEigenVariations::excludeNamedUncertainty error:" << " name " << name << " not allowed" << std::endl;
  } 
  // in case multiple uncertainties should be discarded
  else if (name.back() == '*'){
    std::string temp_name = name.substr(0, name.size()-1); //remove "*"
    std::vector<std::string> uncs = m_cnt->listUncertainties();
    std::vector<std::string> unc_subgroup;
    std::copy_if(uncs.begin(), uncs.end(), back_inserter(unc_subgroup),
         [&temp_name](const std::string& el) {
           return el.compare(0, temp_name.size(), temp_name) == 0;
         });
    if (m_verbose) std::cout <<"Found a group of uncertainties to exclude: " <<name <<" found " <<unc_subgroup.size() <<" uncertainties corresponding to the query" <<std::endl;
    for (const auto& single_name : unc_subgroup){
      // only really add if the entry is not yet in the list
      if (m_namedIndices.find(single_name) == m_namedIndices.end()) {
        if (m_verbose) std::cout << "Name : " << single_name << std::endl;
        m_named.push_back(std::pair<TH1*, TH1*>(0, 0));
        m_namedIndices[single_name] = m_named.size()-1;
      }
    }
  }
  else if (! cnt->GetValue(name.c_str())){
    std::cerr << "CalibrationDataEigenVariations::excludeNamedUncertainty error:" << " uncertainty named " << name << " not found" << std::endl;
  } 
}

excludeNamedUncertainty() [2/2]

void CalibrationDataGlobalEigenVariations::excludeNamedUncertainty	(	const std::string &	name,
		const std::string &	flavour )

Definition at line 1523 of file CalibrationDataEigenVariations.cxx.

{
  // Exclude the source of uncertainty identified by the given name from being used
  // in the construction of the covariance matrix to be diagonalised.
  // Notes:
  // - Some names returned by CalibrationDataContainer::listUncertainties() are not
  //   meaningful in this context, and specifying them is not allowed.
  // - Once the eigenvector diagonalisation has been carried out, this method may
  //   not be used anymore.
 
  if (m_initialized) std::cerr << "CalibrationDataGlobalEigenVariations::excludeNamedUncertainty error:" << " initialization already done" << std::endl;
 
  else if (name == "comment"    || name == "result"   || name == "systematics" || 
           name == "statistics" || name == "combined" || name == "extrapolation" || 
           name == "MCreference" || name == "MChadronisation" || name == "ReducedSets" || 
           name == "excluded_set" || name == "" || name == " ")
    std::cerr << "CalibrationDataGlobalEigenVariations::excludeNamedUncertainty error:" << " name [" << name << "] not allowed" << std::endl;
 
  // Note: we WANT to exclude named uncertainties FROM ALL FLAVOURS, even if they don't contain the uncertainty (just fill with zeros). So I won't do this check here.
  // only really add if the entry is not yet in the list
  else if (m_flav_namedIndices[flavour].find(name) == m_flav_namedIndices[flavour].end()) {
    m_flav_named[flavour].push_back(std::pair<TH1*, TH1*>(0, 0));
    m_flav_namedIndices[flavour][name] = m_flav_named[flavour].size()-1; // store the index that the name variation pair is stored with in m_named
  }
}

getEigenCovarianceMatrix()

TMatrixDSym CalibrationDataGlobalEigenVariations::getEigenCovarianceMatrix ( )

virtual

also provide (some) access to the underlying information: covariance matrix corresponding to eigenvector variations

Reimplemented from Analysis::CalibrationDataEigenVariations.

Definition at line 1385 of file CalibrationDataEigenVariations.cxx.

{
  // Construct the block (global) covariance matrix that is to be diagonalised.
  // Note that extrapolation uncertainties (identified by the keyword "extrapolation,
  // this will pertain mostly to the extrapolation to high jet pt) are always excluded
  // since by definition they will not apply to the normal calibration bins. Instead
  // this uncertainty has to be dealt with as a named variation. In addition there are
  // other items ("combined", "systematics") that will not be dealt with correctly
  // either and hence are excluded as well).
  //
  // Note that if an explicit covariance matrix is supplied (at present this may be
  // the case only for statistical uncertainties: in the case of "continuous tagging",
  // multinomial statistics applies so bin-to-bin correlations exist) this will be
  // used instead of constructing the statistical uncertainties' covariance matrix on
  // the fly.
 
  std::map<std::string, TMatrixDSym> cov_matrices; // temporary store, just to aid in printing out the individual systematic covariance matrices
  TMatrixDSym global_covariance(m_blockmatrixsize);
  // Then loop through the list of (other) uncertainties
  for (const std::string& unc : m_all_shared_systematics){
    std::vector<int> flavs_in_common;
    TString tunc(unc);
    std::vector<double> comb_syst; // this vector combines the TH1 uncertainty bins for each flavour into one object -> stored in comb_systematics for now, but meant for covariance matrix method
    // For the fixed cut case, we want to bin by groups of flavour > pT (i.e. "blocks" of flavour)
    // For the continuous case, we want to bin by groups of flavour > tagweight > pT (i.e. "blocks" of flavour, containing tagweight blocks... more complicated, but works on same principle)
    for (const auto& flavour : m_flavours){
      Analysis::CalibrationDataHistogramContainer* c = m_histcontainers[flavour]; // pointer to the flavour container
      int flavour_size = 0; // Store the length of the uncertainty in the flavour, initialize to zero
      if (c) {
        // Now we want to get the histogram for this systematic, for this flavour
        TH1* hunc = dynamic_cast<TH1*>(c->GetValue(tunc.Data()));
        TH1* ref = dynamic_cast<TH1*>(c->GetValue("result")); // retrieving this just in case the uncertainty doesn't exist for this flavour, just need it to get dimensions right
        if (not ref){
           if (m_verbose) std::cout << " There was no uncertainty OR SF/EFF results... Are you sure you have the right CDIContainer path?" << std::endl;
           continue;
        }
        int tagweightax = c->getTagWeightAxis(); // for handling the continuous case(s)
        if (hunc){
          flavs_in_common.push_back(1); // report that we have this uncertainty in this flavour
          Int_t nbinx = hunc->GetNbinsX(), nbiny = hunc->GetNbinsY(), nbinz = hunc->GetNbinsZ();
          Int_t rows = nbinx;
          if (hunc->GetDimension() > 1) rows *= nbiny;
          if (hunc->GetDimension() > 2) rows *= nbinz;
          flavour_size = rows; // Record the number of bins in the uncertainty TH1
        } else {
          flavs_in_common.push_back(0); // If the uncertainty doesn't exist for that flavour, just report, we'll set to zero in the combined vector
          // Because the uncertainty doesn't exist for this flavour, we just get the dimensions we need
          Int_t nbinx = ref->GetNbinsX(), nbiny = ref->GetNbinsY(), nbinz = ref->GetNbinsZ();
          Int_t rows = nbinx;
          if (ref->GetDimension() > 1) rows *= nbiny;
          if (ref->GetDimension() > 2) rows *= nbinz;
          flavour_size = rows;
        }
        
        // Now we can loop through the bins of the flavour uncertainty, adding them onto the combined systematic
        if (tagweightax == -1){ //<--- i.e. NOT continuous, but is a fixed cut WP
          for (int i = 1 ; i <= flavour_size ; i++){
            if (hunc){
              Int_t bin = hunc->GetBin(1,i,1);
              double unc_val = hunc->GetBinContent(bin);
              comb_syst.push_back(unc_val); // if uncertainty, push uncertainty bin content to combined systematic vector
            } else {
              comb_syst.push_back(0.0); // if no uncertainty, push 0's
            }
          }
        } else if (tagweightax == 0){
          // X axis is the tagweight axis, meaning Y is pt, Z is abs(eta)
          for (Int_t xbins = 1 ; xbins <= ref->GetNbinsX() ; xbins++){
            for (Int_t zbins = 1 ; zbins <= ref->GetNbinsZ() ; zbins++){
              for (Int_t ybins = 1 ; ybins <= ref->GetNbinsY() ; ybins++){
                if (hunc){
                  Int_t bin = hunc->GetBin(xbins,ybins,zbins); 
                  double unc_val = hunc->GetBinContent(bin); 
                  comb_syst.push_back(unc_val); // if uncertainty, push uncertainty bin content to combined systematic vector 
                } else {
                  comb_syst.push_back(0.0); // if no uncertainty, push 0's
                }
                // And now we should be constructing the initial covariance matrix for block continuous correctly
              }
            }
          }
        } else if (tagweightax == 1){
          // Y axis is the tagweight axis, meaning X is pt, Z is abs(eta)
          for (Int_t ybins = 1 ; ybins <= ref->GetNbinsY() ; ybins++){
            for (Int_t zbins = 1 ; zbins <= ref->GetNbinsZ() ; zbins++){
              for (Int_t xbins = 1 ; xbins <= ref->GetNbinsX() ; xbins++){
                if (hunc){
                  Int_t bin = hunc->GetBin(xbins,ybins,zbins); 
                  double unc_val = hunc->GetBinContent(bin);
                  comb_syst.push_back(unc_val); // if uncertainty, push uncertainty bin content to combined systematic vector 
                } else {
                  comb_syst.push_back(0.0); // if no uncertainty, push 0's
                }
                // And now we should be constructing the initial covariance matrix for block continuous correctly
              }
            }
          }
        } else if (tagweightax == 2){
          // Z axis is the tagweight axis, meaning X is pt, Y is abs(eta)
          for (Int_t zbins = 1 ; zbins <= ref->GetNbinsZ() ; zbins++){
            for (Int_t ybins = 1 ; ybins <= ref->GetNbinsY() ; ybins++){
              for (Int_t xbins = 1 ; xbins <= ref->GetNbinsX() ; xbins++){
                if (hunc){
                  Int_t bin = hunc->GetBin(xbins,ybins,zbins); 
                  double unc_val = hunc->GetBinContent(bin); 
                  comb_syst.push_back(unc_val); // if uncertainty, push uncertainty bin content to combined systematic vector 
                } else {
                  comb_syst.push_back(0.0); // if no uncertainty, push 0's
                }
                // And now we should be constructing the initial covariance matrix for block continuous correctly
              }
            }
          }
        }
      }
    }
    //comb_systematics.insert({unc, comb_syst}); // after having looped through the bins for each flavour, the comb_syst vector should be completed and added (or rather, made into a covariance matrix)
    TMatrixDSym unc_cov(comb_syst.size());
    if (unc == "statistics"){
      unc_cov = getStatCovarianceMatrix(comb_syst, m_statVariations); // we want to handle the "statistics" uncertainty differently
    } else {
      unc_cov = getSystCovarianceMatrix(comb_syst);
    }
    cov_matrices.insert({unc, unc_cov}); // add the covariance matrix for the uncertainty to this map, this is temporary for testing/plotting purposes
    
    // To look only at uncertainties that pertain to more than one flavour (2, 3 or all 4) we can store the covariances separately for later saving and plotting
    if (flavs_in_common[0]+flavs_in_common[1]+flavs_in_common[2]+flavs_in_common[3] > 1){
      m_only_shared_systematics.insert(unc) ; // mark the shared systematics...
    }
 
    global_covariance += unc_cov;
  }
  return global_covariance;
}

getEigenCovarianceMatrixFromVariations()

TMatrixDSym CalibrationDataEigenVariations::getEigenCovarianceMatrixFromVariations ( )

inherited

covariance matrix corresponding to eigenvector variations constructed from the eigen-variation

Definition at line 473 of file CalibrationDataEigenVariations.cxx.

{
  // Construct the (Eigen-variation part of the) covariance matrix from the individual variations.
  // This must be called _after_ initialize()!
 
  // retrieve the central calibration
  TH1        *result = dynamic_cast<TH1*>(m_cnt->GetValue("result")); 
  if (!result){
    std::cerr<<"CalibrationDataEigenVariations::getEigenCovarianceMatrixFromVariations(): dynamic cast failed\n";
    return TMatrixDSym();
  }
  TMatrixD    jac = getJacobianReductionMatrix();
  int         nbins = jac.GetNcols();
  TMatrixDSym cov(nbins);
  auto variation = std::make_unique<double[]>(nbins);
 
  for (const std::pair<TH1*, TH1*>& it : m_eigen){ // m_eigen is vector of pairs of TH1* which point to TH1's representing the upVariation and downVariation respectively
    TH1* resultVariedUp = it.first; // <--------------- This is the "result" but "varied" upwards - i.e. how the result would look like if the systematic "it" was imposed
    for (unsigned int u = 0; u < (unsigned int) nbins; ++u) variation[u] = resultVariedUp->GetBinContent(u) - result->GetBinContent(u);
    for (int u = 0; u < nbins; ++u){ 
      for (int v = 0; v < nbins; ++v){
        cov(u, v) += variation[u]*variation[v];
      }
    }
  }
 
  return cov;
}

getEigenvectorVariation() [1/2]

bool CalibrationDataEigenVariations::getEigenvectorVariation ( unsigned int variation, TH1 *& up, TH1 *& down )

inherited

obtain the "up" and "down" variations for the given eigenvector number.

The return value will be false if the eigenvector number is invalid.

Definition at line 1022 of file CalibrationDataEigenVariations.cxx.

{
  // Return the pointers to the up- and downward variation histogram for the specified
  // eigenvector variation. In case of an invalid variation number, null pointers will
  // be returned and the return value will be false.
  //
  //     variation:   eigenvector variation number
  //     up:          (reference to) pointer to upward variation histogram
  //     down:        (reference to) pointer to downward variation histogram
 
  if (! m_initialized) initialize();
  
  if (variation < m_eigen.size()) {
    up   = m_eigen[variation].first;
    down = m_eigen[variation].second;
    return true;
  }
 
  up = down = 0;
  return false;
}

getEigenvectorVariation() [2/2]

bool CalibrationDataGlobalEigenVariations::getEigenvectorVariation	(	const std::string &	flavour,
		unsigned int	variation,
		TH1 *&	up,
		TH1 *&	down )

Definition at line 1961 of file CalibrationDataEigenVariations.cxx.

{
  // Return the pointers to the up- and downward variation histogram for the specified
  // eigenvector variation. In case of an invalid variation number, null pointers will
  // be returned and the return value will be false.
  //
  //     variation:   eigenvector variation number
  //     up:          (reference to) pointer to upward variation histogram
  //     down:        (reference to) pointer to downward variation histogram
 
  if (! m_initialized) initialize();
  std::vector<std::pair<TH1*,TH1*>> flav_variations = m_flav_eigen.find(flavour)->second;
  if (variation < flav_variations.size()){
    up   = flav_variations[variation].first;
    down = flav_variations[variation].second;
    return true;
  }
 
  up = down = 0;
  return false;
}

getJacobianReductionMatrix() [1/2]

TMatrixD CalibrationDataEigenVariations::getJacobianReductionMatrix ( )

inherited

matrix to remove unecessary rows and columns from covariance

Definition at line 504 of file CalibrationDataEigenVariations.cxx.

{
  // Construct the matrix that removes the rows and columns that fail to influence
  // the eigen-variations. To reduce the covariance matrix, do the following:
  // 
  // TMatrixDSym cov = getEigenCovarianceMatrix();
  // TMatrixDSym jac = getJacobianReductionMatrix(); // <------------ This is an upper triangular matrix of 1's. Similarity transformation will do...
  // TMatrixDSym redSystematicCovMatrix = cov.Similarity(jac);
 
  // retrieve the central calibration
  TH1* result = dynamic_cast<TH1*>(m_cnt->GetValue("result"));
  if (not result){
    std::cerr<<"CalibrationDataEigenVariations::getJacobianReductionMatrix(): dynamic cast failed\n";
    return TMatrixD();
  }
 
  // loop over the uncertainties to construct the covariance matrix for all uncertainties
  // to be addressed using the eigenvector method.
 
  // Retrieve the un-compressed Eigenvector variation covariance matrix
  // (only needed to check for unexpected singularities)
  TMatrixDSym cov = getEigenCovarianceMatrix();
 
  // Compute the original number of bins
  int nbins = result->GetNbinsX()+2;
  int ndim  = result->GetDimension();
  if (ndim > 1) nbins*= (result->GetNbinsY()+2);
  if (ndim > 2) nbins*= (result->GetNbinsZ()+2);
 
  // Start by "compressing" the covariance matrix (removing columns/rows containing zeros only)
  int nZeros=0;
  std::vector<int> zeroComponents;
  if (cov.GetNrows() != nbins) {
    std::cerr << " error: covariance matrix size (" << cov.GetNrows() << ") doesn't match histogram size (" << nbins << ")" << std::endl;
    return TMatrixDSym();
  }
 
  // First flag all the zeros
  for (int i = 0; i < nbins; ++i) {
    // Directly identify the under- and overflow bins
    Int_t binx, biny, binz;
    result->GetBinXYZ(i, binx, biny, binz);
    if ((binx == 0 || binx == result->GetNbinsX()+1) || 
       (ndim > 1 && (biny == 0 || biny == result->GetNbinsY()+1)) || 
       (ndim > 2 && (binz == 0 || binz == result->GetNbinsZ()+1))) 
    {
      ++nZeros;
      zeroComponents.push_back(i);
    }
    // Try a first (quick) identification of rows/columns of zeros by the first element in each row
    // If "successful", check the whole row in more detail
    else if (fabs(cov(i,0)) < 1e-10) {
      bool isThereANonZero(false);
      for (int j = 0; j < nbins; ++j) {
        if (fabs(cov(i,j)) > 1e-10) {
          isThereANonZero=true; break;
        }
      }
      if (! isThereANonZero) {
        ++nZeros;
        zeroComponents.push_back(i);
      }
    }
  }
 
  // **** COMMENTED OUT FOR NOW. 
  // Leave it here in case the calibration method will change again in the future. 
  // No need to reweight the SF by the efficiency of that bin (MCreference always = 0)
 
  // Determine whether the container is for "continuous" calibration.
  // This is important since the number of independent scale factors (for each pt or eta bin)
  // is reduced by 1 compared to the number of tag weight bins (related to the fact that the fractions
  // of events in tag weight bins have to sum up to unity).
  // int axis = m_cnt->getTagWeightAxis();
  // bool doContinuous = false;  unsigned int weightAxis = 0;
 
  // if (axis >= 0) {
  //   doContinuous = true;
  //   // weightAxis = (unsigned int) axis;
  //   // In this case, verify that the special "uncertainty" entry that is in fact the reference MC tag
  //   // weight fractions is present. These tag weight fractions are needed in order to carry out the
  //   // diagonalisation successfully.
  //   if (! dynamic_cast<TH1*>(m_cnt->GetValue("MCreference"))) {
  //     std::cerr << " Problem: continuous calibration object found without MC reference tag weight histogram " << std::endl;
  //     return TMatrixDSym();
  //   }
  // }
 
  // Only relevant for continuous calibration containers, but in order to void re-computation we
  // retrieve them here
  // Int_t nbinx = result->GetNbinsX()+2, nbiny = result->GetNbinsY()+2, nbinz = result->GetNbinsZ()+2;
 
  // // If we are indeed dealing with a "continuous" calibration container, ignore one tag weight row
  // const int skipTagWeightBin = 1; // NB this follows the histogram's bin numbering
  // if (doContinuous) {
  //   for (Int_t binx = 1; binx < nbinx-1; ++binx)
  //     for (Int_t biny = 1; biny < nbiny-1; ++biny)
  //       for (Int_t binz = 1; binz < nbinz-1; ++binz) {
  //         if ((weightAxis == 0 && binx == skipTagWeightBin) ||
  //             (weightAxis == 1 && biny == skipTagWeightBin) ||
  //             (weightAxis == 2 && binz == skipTagWeightBin)) {
  //           // At this point we simply add these to the 'null' elements
  //           ++nZeros;
  //           zeroComponents.push_back(result->GetBin(binx, biny, binz));
  //         }
  //       }
  // }
 
  if (nZeros >= nbins) {
    std::cerr << " Problem. Found n. " << nZeros << " while size of matrix is " << nbins << std::endl;
    return TMatrixDSym();
  }
 
  int size = nbins - nZeros; // <--- Tis the size of the matrix removing zeros from the covariance matrix
 
  TMatrixT<double> matrixVariationJacobian(size,nbins);
  int nMissed=0;
  for (int i = 0; i < nbins; ++i) { //full size
    bool missed=false;
    for (unsigned int s = 0; s < zeroComponents.size(); ++s) { // <-------- Basically what this does is it flags "missed" for a given "i" of the full bin size
      if (zeroComponents.at(s) == i) {                        //  <-------- if "i" is in "zeroComponents". Breaks (because it found that it's to be missed)
        missed = true;
        break;
      }
    }
    if (missed) {                                       //        <-------- Finally, if "i" is to be missed, increase "nMissed" by one, and....
      ++nMissed;
      continue; 
    }
    matrixVariationJacobian(i-nMissed,i)=1;         //            <-------- ... this ALWAYS adds a one. If zero "nMissed", add to diagonal. otherwise off-diagonal
  }                                                //             <-------- This matrix only add 1's on/off diagonal, upper triangular matrix
 
  return matrixVariationJacobian;
}

getJacobianReductionMatrix() [2/2]

TMatrixD CalibrationDataGlobalEigenVariations::getJacobianReductionMatrix

(

TMatrixDSym &

cov

)

Definition at line 1883 of file CalibrationDataEigenVariations.cxx.

{
  // Construct the matrix that removes the rows and columns that fail to influence
  // the eigen-variations. To reduce the covariance matrix, do the following:
  // Note: Previous version called "getEigenCovariance" whereas here we pass it in by reference to save on compute
  // This "cov" is the "uncompressed" eigenvariation covariance matrix for all uncertainties.
  // We will subsequently compress it (removing columns/rows containing zeros only) 
  // and then construct the hopefully rectangular matrix that can remove these ones from the covariance matrix
 
  int nZeros = 0;
  std::vector<int> zeroComponents;
 
  // First flag all the zeros
  for (int i = 0 ; i < m_blockmatrixsize ; i++){
    if (fabs(cov(i,0)) < 1e-10){
      bool isThereANonZero = false;
      for (int j = 0 ; j < m_blockmatrixsize ; j++){
        // Try a first (quick) identification of rows/columns of zeros by the first element in each row
        // If "successful", check the whole row in more detail
        if (fabs(cov(i,j)) > 1e-10){
          isThereANonZero = true;
          break;
        }
      }
      if (! isThereANonZero){
        ++nZeros;
        zeroComponents.push_back(i) ;
      }
    }
  }
 
  if (nZeros >= m_blockmatrixsize) {
    std::cerr << " Problem. Found n. " << nZeros << " while size of matrix is " << m_blockmatrixsize << std::endl;
    return TMatrixDSym();
  }
  
  int size = m_blockmatrixsize - nZeros;
  TMatrixT<double> matrixVariationJacobian(size,m_blockmatrixsize);
  int nMissed = 0;
 
  for (int i = 0; i < m_blockmatrixsize; ++i) { // full size
    bool missed = false;
    for (unsigned int s = 0 ; s < zeroComponents.size(); ++s) { // <--- Basically what this does is it flags "missed" for a given "i" of the full bin size
      if (zeroComponents.at(s) == i) { // <--- if "i" is in "zeroComponents". Breaks (because it found that it's to be missed)
        missed = true;
        break;
      }
    }
    if (missed) {  // <-------- Finally, if "i" is to be missed, increase "nMissed" by one, and....
      ++nMissed;
      continue;
    }
 
    matrixVariationJacobian(i-nMissed,i)=1; //  <-------- ... this ALWAYS adds a one. If zero "nMissed", add to diagonal. otherwise off-diagonal
  }
 
  return matrixVariationJacobian;
}

getNamedVariation() [1/4]

bool CalibrationDataEigenVariations::getNamedVariation ( const std::string & name, TH1 *& up, TH1 *& down )

inherited

obtain the "up" and "down" variations for the named uncertainty.

The return value will be false if the given name is not listed as being excluded from the eigenvector calculations.

Definition at line 1047 of file CalibrationDataEigenVariations.cxx.

{
  // Return the pointers to the up- and downward variation histogram for the specified
  // named variation. In case of an invalid named variation, null pointers will
  // be returned and the return value will be false.
  //
  //     name:        named variation
  //     up:          (reference to) pointer to upward variation histogram
  //     down:        (reference to) pointer to downward variation histogram
 
  map<string, unsigned int>::const_iterator it = m_namedIndices.find(name);
  if (it != m_namedIndices.end()) return getNamedVariation(it->second, up, down);
 
  up = down = 0;
  return false;
}

getNamedVariation() [2/4]

bool CalibrationDataEigenVariations::getNamedVariation ( unsigned int nameIndex, TH1 *& up, TH1 *& down )

inherited

obtain the "up" and "down" variations for the source uncertainty pointed to by the given index (which is assumed to correspond to the value retrieved using getNamedVariationIndex()).

The return value will be false if the index is out of bounds.

Definition at line 1067 of file CalibrationDataEigenVariations.cxx.

{
  // Return the pointers to the up- and downward variation histogram for the specified
  // named variation. In case of an invalid named variation number, null pointers will
  // be returned and the return value will be false.
  //
  //     nameIndex:   named variation number
  //     up:          (reference to) pointer to upward variation histogram
  //     down:        (reference to) pointer to downward variation histogram
 
  if (! m_initialized) initialize();
 
  if (nameIndex < m_named.size()) {
    up   = m_named[nameIndex].first;
    down = m_named[nameIndex].second;
    return true;
  }
 
  up = down = 0;
  return false;
}

getNamedVariation() [3/4]

bool CalibrationDataGlobalEigenVariations::getNamedVariation	(	const std::string &	flavour,
		const std::string &	name,
		TH1 *&	up,
		TH1 *&	down )

Definition at line 1985 of file CalibrationDataEigenVariations.cxx.

{
  // Return the pointers to the up- and downward variation histogram for the specified
  // named variation. In case of an invalid named variation, null pointers will
  // be returned and the return value will be false.
  //
  //     name:        named variation
  //     up:          (reference to) pointer to upward variation histogram
  //     down:        (reference to) pointer to downward variation histogram
 
  // Find the named variation index (if it exists) and pass to "getNamedVariation"
  std::map<std::string, unsigned int>::const_iterator it = (m_flav_namedIndices.find(flavour)->second).find(name);
  if (it != (m_flav_namedIndices.find(flavour)->second).end()) return getNamedVariation(it->second, flavour, up, down);
 
  up = down = 0;
  return false;
}

getNamedVariation() [4/4]

bool CalibrationDataGlobalEigenVariations::getNamedVariation	(	unsigned int	nameIndex,
		const std::string &	flavour,
		TH1 *&	up,
		TH1 *&	down )

Definition at line 2005 of file CalibrationDataEigenVariations.cxx.

{
  // Return the pointers to the up- and downward variation histogram for the specified
  // named variation. In case of an invalid named variation number, null pointers will
  // be returned and the return value will be false.
  //
  //     nameIndex:   named variation number
  //     up:          (reference to) pointer to upward variation histogram
  //     down:        (reference to) pointer to downward variation histogram
 
  if (! m_initialized) initialize();
 
  if (nameIndex < m_flav_named[flavour].size()) {
    up   = m_flav_named[flavour][nameIndex].first;
    down = m_flav_named[flavour][nameIndex].second;
    return true;
  }
 
  up = down = 0;
  return false;
}

getNamedVariationIndex() [1/2]

unsigned int CalibrationDataEigenVariations::getNamedVariationIndex ( const std::string & name ) const

inherited

retrieve the integer index corresponding to the named variation.

This can be used to speed up computations by avoiding string comparisons. Note that this function is not protected against passing an invalid name.

Definition at line 1092 of file CalibrationDataEigenVariations.cxx.

{
  // Return the integer index corresponding to the named variation.
  // Note that no checks are made on the validity of the name provided.
 
  map<string, unsigned int>::const_iterator it = m_namedIndices.find(name);
  return it->second;
}

getNamedVariationIndex() [2/2]

unsigned int CalibrationDataGlobalEigenVariations::getNamedVariationIndex	(	const std::string &	name,
		const std::string &	flavour ) const

Definition at line 2030 of file CalibrationDataEigenVariations.cxx.

{
  // Return the integer index corresponding to the named variation.
  // Note that no checks are made on the validity of the name provided.
  std::map<std::string, unsigned int>::const_iterator it = (m_flav_namedIndices.find(flavour)->second).find(name);
  return it->second;
}

getNumberOfEigenVariations() [1/2]

unsigned int CalibrationDataEigenVariations::getNumberOfEigenVariations ( )

inherited

retrieve the number of eigenvector variations

Definition at line 1014 of file CalibrationDataEigenVariations.cxx.

{
  if (! m_initialized) initialize();
  return m_eigen.size();
}

getNumberOfEigenVariations() [2/2]

unsigned int CalibrationDataGlobalEigenVariations::getNumberOfEigenVariations

(

const std::string &

flavour

)

Definition at line 1948 of file CalibrationDataEigenVariations.cxx.

{
  if (! m_initialized) initialize();
  
  if (m_flav_eigen.find(flavour) == m_flav_eigen.end()){
    if (m_verbose) std::cout << "No m_flav_eigen for flavour " << flavour << std::endl;
    return 0;
  }
  return (m_flav_eigen.find(flavour)->second).size();
}

getNumberOfNamedVariations()

unsigned int CalibrationDataEigenVariations::getNumberOfNamedVariations ( ) const

inherited

retrieve the number of named variations

Definition at line 992 of file CalibrationDataEigenVariations.cxx.

{
  // Returns the number of named variations
 
  return m_namedIndices.size();
}

initialize()

void CalibrationDataGlobalEigenVariations::initialize ( double min_variance = 1.0E-6 )

virtual

carry out the eigenvector computations.

Exclusion of any source of uncertainty has to be done before calling this method

Reimplemented from Analysis::CalibrationDataEigenVariations.

Definition at line 1566 of file CalibrationDataEigenVariations.cxx.

{
  // This is this class's most important method, in the sense that it does all the
  // math and constructs all "variation" histograms (for both eigenvector and named
  // named variations). This constitutes the full initialisation of the class.
  // This method is meant to be called only after all calls (if any) to the
  // CalibrationDataGlobalEigenVariations::excludeNamedUncertainty() method.
 
 
  // First step: construct the block covariance matrix
  TMatrixDSym cov = getEigenCovarianceMatrix(); 
 
  TMatrixDSym corr(cov); // We want to construct the correlation matrix in order to compare the final eigenvariations correlation matrix to it
  for (int row = 0 ; row < cov.GetNrows() ; row++){
    for (int col = 0 ; col < cov.GetNcols() ; col++){
      double rowvar = sqrt(cov(row,row));
      double colvar = sqrt(cov(col,col));
      corr(row,col) = corr(row,col)/(rowvar * colvar); // divide each element by the variance
    }
  }

 
  // Second step: create the variations for the named sources of uncertainty (easy...)
  // News flash: This isn't so easy now, as it has to be done for all the flavours to which the uncertainty pertains
  // the following are temporary data structures
  std::vector<double> combined_result; 
  std::map<std::string, int> flav_bins; // store the nbins of each flavour histogram for later use in "reporting"
  // and eventually, at the end, we will want to separate the flavour blocks out...
  for (const auto& flavour : m_flavours) {
    // for each flavour, we want to combine the results and uncertainty values across flavours into one "block" vector
    // retrieve the central calibration
    Analysis::CalibrationDataHistogramContainer* c = m_histcontainers[flavour]; // pointer to the flavour container
 
    TH1* result = dynamic_cast<TH1*>(c->GetValue("result"));
    if (not result){
      std::cerr<<"CalibrationDataGlobalEigenVariations::initialize: dynamic cast failed\n ";
      continue;
    }
    //construct the combined_result and combined_named_variations (up and down)
    if (c->getTagWeightAxis() == -1){ // For fixed cut WP, the Y axis **should** be the pT axis (but can it can potentially be different in the future)
      flav_bins[flavour] = result->GetNbinsY(); // Add the number of bins of the result histogram with non-zero results...  
      if (m_verbose) std::cout << "flav_bins["<<flavour<<"] = " << flav_bins[flavour] << std::endl;
      for(int i = 0 ; i < flav_bins[flavour] ; i++){
        // push bin content onto the combined_result vector
        Int_t bin = result->GetBin(1,i+1); // This will only pick up the CONTENTS, not of the underflow bin
        double res_value = result->GetBinContent(bin);
        combined_result.push_back(res_value);
      }
    } else if (c->getTagWeightAxis() == 0) { // for continuous WP, the taxweight axis determines which axis is pt and |eta|
      flav_bins[flavour] = result->GetNbinsX()*result->GetNbinsY();
      int tagbins = result->GetNbinsX();
      int ptbins = result->GetNbinsY();
      if (m_verbose) std::cout << "flav_bins["<<flavour<<"] = " << flav_bins[flavour] << std::endl;
      for(int i = 0 ; i < tagbins ; i++){
        for(int j = 0 ; j < ptbins ; j++){
          // push bin content onto the combined_result vector
          Int_t bin = result->GetBin(j+1,1,i+1); // This will only pick up the CONTENTS, not of the underflow bin
          double res_value = result->GetBinContent(bin);
          combined_result.push_back(res_value);
        }
      }
    } else if (c->getTagWeightAxis() == 1) {
      flav_bins[flavour] = result->GetNbinsX()*result->GetNbinsY();
      int tagbins = result->GetNbinsY();
      int ptbins = result->GetNbinsX();
      if (m_verbose) std::cout << "flav_bins["<<flavour<<"] = " << flav_bins[flavour] << std::endl;
      for(int i = 0 ; i < tagbins ; i++){
        for(int j = 0 ; j < ptbins ; j++){
          // push bin content onto the combined_result vector
          Int_t bin = result->GetBin(j+1,1,i+1); // This will only pick up the CONTENTS, not of the underflow bin
          double res_value = result->GetBinContent(bin);
          combined_result.push_back(res_value);
        }
      }
    } else if (c->getTagWeightAxis() == 2) {
      flav_bins[flavour] = result->GetNbinsX()*result->GetNbinsZ();
      int tagbins = result->GetNbinsZ();
      int ptbins = result->GetNbinsX();
      if (m_verbose) std::cout << "flav_bins["<<flavour<<"] = " << flav_bins[flavour] << std::endl;
      for(int i = 0 ; i < tagbins ; i++){
        for(int j = 0 ; j < ptbins ; j++){
          // push bin content onto the combined_result vector
          Int_t bin = result->GetBin(j+1,1,i+1); // This will only pick up the CONTENTS, not of the underflow bin
          double res_value = result->GetBinContent(bin);
          combined_result.push_back(res_value);
        }
      }
    }
 
 
    // loop over the excluded uncertainties for this flavour, and construct their actual variations...
    // the "m_flav_namedIndices" are constructed within "excludeNamedUncertainties", which is called in the constructor
    for (map<string, unsigned int>::iterator it = m_flav_namedIndices[flavour].begin(); it != m_flav_namedIndices[flavour].end(); ++it) {
      TH1* hunc = (TH1*) c->GetValue(it->first.c_str()); // this should store the name uncertainty, if it exists for this flavour 
      // I need to test if this uncertainty actually exists, and if it doesn't just use a ZERO variation histogram explicitly
      // but even if it doesn't exist, we want to have it, so that the indices match between flavours
      pair<TH1*, TH1*>& p = m_flav_named[flavour][it->second];
      TString namedvar("namedVar");
      namedvar += it->first.c_str();
      TString namedvarUp(namedvar);   namedvarUp   += "_up";
      TString namedvarDown(namedvar); namedvarDown += "_down";
      TH1* resultVariedUp   = (TH1*)result->Clone(namedvarUp);   
      TH1* resultVariedDown = (TH1*)result->Clone(namedvarDown); 
      if (hunc){
        resultVariedUp->Add(hunc, 1.0);    resultVariedUp->SetDirectory(0);
        resultVariedDown->Add(hunc, -1.0); resultVariedDown->SetDirectory(0);
      } else {
        resultVariedUp->SetDirectory(0);
        resultVariedDown->SetDirectory(0);
      }
      p.first  = resultVariedUp;
      p.second = resultVariedDown;
    } // End of uncertainty in flavour loop 
    
    // Now handle the extrapolation uncertainties per flavour...
    // Refinement: add the "extrapolation" uncertainty as a named uncertainty, if the histogram is provided
    // This is a bit special, since the extrapolation uncertainty histogram has a different size than other histograms
    if (TH1* hunc = (TH1*)c->GetValue("extrapolation")) { // this is just saying "if it exists"...
      TH1* resultVariedUp   = (TH1*) hunc->Clone("extrapolation_up");   resultVariedUp->SetDirectory(0);
      TH1* resultVariedDown = (TH1*) hunc->Clone("extrapolation_down"); resultVariedDown->SetDirectory(0);
      Int_t nbinx = hunc->GetNbinsX()+2, nbiny = hunc->GetNbinsY()+2, nbinz = hunc->GetNbinsZ()+2;
      Int_t firstbinx = hunc->GetXaxis()->FindFixBin(result->GetXaxis()->GetBinCenter(1));
      Int_t firstbiny = result->GetDimension() > 1 ? hunc->GetYaxis()->FindFixBin(result->GetYaxis()->GetBinCenter(1)) : 1;
      Int_t firstbinz = result->GetDimension() > 2 ? hunc->GetZaxis()->FindFixBin(result->GetZaxis()->GetBinCenter(1)) : 1;
      for (Int_t binx = 1; binx < nbinx-1; ++binx) {
        Int_t binxResult = binx - firstbinx + 1;
        if (binxResult < 1) binxResult = 1;
        if (binxResult > result->GetNbinsX()) binxResult = result->GetNbinsX();
        for (Int_t biny = 1; biny < nbiny-1; ++biny) {
          Int_t binyResult = biny - firstbiny + 1;
          if (binyResult > result->GetNbinsY()) binyResult = result->GetNbinsY();
          for (Int_t binz = 1; binz < nbinz-1; ++binz) {
            Int_t binzResult = binz - firstbinz + 1;
            if (binzResult > result->GetNbinsZ()) binzResult = result->GetNbinsZ();
            Int_t bin = hunc->GetBin(binx, biny, binz);
            double contResult = result->GetBinContent(binxResult, binyResult, binzResult);
            resultVariedUp->SetBinContent(bin, contResult + hunc->GetBinContent(bin));
            resultVariedDown->SetBinContent(bin, contResult - hunc->GetBinError(bin));
          }
        }
      }
      m_flav_named[flavour].push_back(std::make_pair(resultVariedUp, resultVariedDown));
      m_flav_namedExtrapolation[flavour] = m_flav_namedIndices[flavour]["extrapolation"] = m_flav_named[flavour].size()-1;
    }
 
 
  } // End flavour loop  
  
 
 
  // Third step: compute the eigenvector variations corresponding to the remaining sources of uncertainty
  // First, build the combined_result vector into a TH1
  std::unique_ptr<TH1> comb_result(new TH1D("combined_result", "", combined_result.size(), 0., 1.));
  int nbins = comb_result->GetNbinsX()+2;
  int ndim  = comb_result->GetDimension();
  if (ndim > 1) nbins*= (comb_result->GetNbinsY()+2);
  if (ndim > 2) nbins*= (comb_result->GetNbinsZ()+2);
 
  // I want to take the contents of "combined_result" and fill in "comb_result" without the under/overflow
  for (unsigned int i=0 ; i<combined_result.size() ; i++){
    // assuming dimension == 1, which should be the case...
    comb_result->SetBinContent(i+1, combined_result[i]); // setting i+1 so we don't start with the underflow bin
  } 
 
  // Get the portions of the covariance matrix that aren't zero, this next step 
  TMatrixT<double> matrixVariationJacobian = getJacobianReductionMatrix(cov); // we pass the covariance matrix in as a starting point
 
  int size = matrixVariationJacobian.GetNrows();
 
  // Reduce the matrix to one without the zeros, using a "similarity" transformation
  const TMatrixDSym matrixCovariance = cov.Similarity(matrixVariationJacobian); // <--- This step removes the zeros
 
  // Carry out the Eigenvector decomposition on this matrix
  TMatrixDSymEigen eigenValueMaker (matrixCovariance);
  TVectorT<double> eigenValues   = eigenValueMaker.GetEigenValues(); 
  TMatrixT<double> eigenVectors  = eigenValueMaker.GetEigenVectors(); 
  TMatrixT<double> matrixVariations (size,size);                     
 
  //compute the total variance by summing the eigenvalues
  m_totalvariance = eigenValues.Sum();
 
  for (int i = 0; i < size; ++i) {
    for (int r = 0; r < size; ++r) {
      //first index is the variation number, second corresponds to the pT bin      // The "eigenvariations" matrix is the "C" matrix which has CKC^T = I with "K" being the o.g. covariance matrix, and "I" is the identity.
      matrixVariations(i,r) = -1.0*eigenVectors[r][i]*sqrt(fabs(eigenValues[i])); //  So the result is a matrix (eigenvariation) which is the eigenvector scaled by the sqrt(eigenvalue)
    } 
  } // <------- matrixVariations: each row is one variation, each column is the pT bin.
 
  TMatrixT<double> matrixVariationsWithZeros = matrixVariations * matrixVariationJacobian; // This step adds in the zero rows again
  
  // Construct the initial set of variations from this
  for (int i = 0; i < matrixVariationsWithZeros.GetNrows(); ++i) {
    TString superstring("eigenVar");
    superstring+=i; 
 
    TString nameUp(superstring);   nameUp   += "_up"; // In the end you get something like "eigenVar5_up"
    TString nameDown(superstring); nameDown += "_down";
    // TString nameUnc(superstring);  nameUnc+= "_unc";
 
    TH1* resultVariedUp   = (TH1*)comb_result->Clone(nameUp);   resultVariedUp->SetDirectory(0);
    TH1* resultVariedDown = (TH1*)comb_result->Clone(nameDown); resultVariedDown->SetDirectory(0);
 
    for (int u = 0; u < comb_result->GetNbinsX(); ++u) {
      resultVariedUp->SetBinContent(u,(comb_result->GetBinContent(u) + matrixVariationsWithZeros(i,u)));
      resultVariedDown->SetBinContent(u,(comb_result->GetBinContent(u) - matrixVariationsWithZeros(i,u)));
    }
 
    m_eigen.push_back(std::make_pair(resultVariedUp, resultVariedDown)); //<--- This is currently storing the FULL/combined variations, which aren't binned with proper bin widths etc.
    // The "proper binning" isn't necessary for pruning purposes. Later on, after pruning, separate the flavour blocks of each eigenvariation and construct the variations for each flavour, storing results in m_flav_eigen
 
  } //end eigenvector size
 
  
  // Without any changes, all eigenvariations are kept due to the variation being non-negligible for SOME flavour...
  // Step 4 : Time for PRUNING
  //  > Check the variation bins to see if any exceed a given threshold value of "min_variance"
  //  > Value of E-6 seems to work well for SFGlobalEigen (found through a scan of possible thresholds)
  //    but this is not necessarily going to hold for future CDI's. This needs to be checked through the
  //    "validate_reduction" function, which can be used to compare how well the reduction captures the total covariance.
  
  // Remove variations that are below the given tolerance (effectively meaning that they don't have any effect)
  IndexSet final_set;
  size_t current_set = 0;
 
  // We set the custom min_variance here 
  min_variance = 1.0E-6;
  for (size_t index = 0; index < m_eigen.size(); ++index) {
    bool keep_variation = false; // guilty until proven innocent
    TH1* up = static_cast<TH1*>(m_eigen[index].first->Clone()); up->SetDirectory(0); // clone the up variation and check it out
    up->Add(comb_result.get(), -1.0); // now we're left with decimal values centered around 0, i.e. 0.02 or -0.02
 
    for (int bin = 1; bin <= nbins; ++bin) {
      if (fabs(up->GetBinContent(bin)) > min_variance) { // If you find even ONE bin with big enough variance, we keep the whole systematic.
        keep_variation = true;
        break;
      }
    }
    if (!keep_variation){ // At this stage, if we find no bins in the systematic with large enough variation, we insert it to "final_set" for removal/pruning
      final_set.insert(current_set);
    } else {
      m_capturedvariance += eigenValues[index];
    }
    delete up;
    ++current_set;
  }
  if (final_set.size() > 0){ // at this point, a set of the indices of negligible variations should have been gathered, proceed to remove them...
    if (m_verbose) std::cout << "CalibrationDataEigenVariations: Removing " << final_set.size() << " eigenvector variations leading to sub-tolerance effects, retaining " << m_eigen.size()-final_set.size() << " variations" << std::endl;
  }
  
  CalibrationDataEigenVariations::removeVariations(final_set); // This method actually performs the reduction. The above logic simply flags which variations to get rid of, inserting them into "final_set"
  
  // AT THIS STAGE: Pruning has already occurred, leaving us with a set of combined eigenvariations in "m_eigen", which we can thence recombine and compute the correlation matrix
  // That correlation matrix can then be compared to the original correlation matrix "corr", simply subtracting one from the other. Better approximations will have close to zero deviation.
  // What follows is the construction of this comparison, and the reporting of the comparison (saving it to file)...  
  std::streamsize ss = std::cout.precision();
  if (m_verbose) std::cout << " The total variance is " << m_totalvariance << " and the reduction captured " << std::setprecision(9) << 100.0*(m_capturedvariance/m_totalvariance) << "% of this." << std::endl;
  std::cout.precision(ss); //restore precision
  if (m_validate){ // <---- this flag can be set manually in a local checkout of this package (temporary fix)
    validate_reduction(m_blockmatrixsize, corr, m_eigen, comb_result.get(), "SFGlobalEigen", m_taggername, m_wp, m_jetauthor); // This method is simply here to report the matrix comparison, for reduction scheme validation
  }
 

  // Let's separate out the flavour histograms now, this should be far simpler to ensure the binning is correct for each flavour.
  // At this stage, we should have constructed a set of combined eigenvariations stored in m_eigen. 
  // Below, we take them one by one, separate, and store the flavour variations in m_flav_eigen, with the proper flavour heading
  // <----- The "mergeVariations" code originally merged the variations in m_eigen. But now that wouldn't work, since m_eigen stores
  //        the combined histograms. We should instead merge the "reported" flavour histograms, as they're properly binned (and physically meaningful)
  //        So if we construct the flavour variations first, keeping the indices all the same, then merge them.
  
 
  for(const std::pair<TH1*,TH1*>& var : m_eigen){
    // now we make use of the same old flavour loop and impose the flavour variation to the flavour container
    TString eigenvarup = var.first->GetName();
    TString eigenvardown = var.second->GetName();
    int bin_baseline = 0; // increment this by flav_bins after getting each flavour block
    for (const std::string& flavour : m_flavours){
      Analysis::CalibrationDataHistogramContainer* c = m_histcontainers[flavour];
      TH1* result = dynamic_cast<TH1*>(c->GetValue("result"));
      if (not result){
        std::cerr<<"CalibrationDataGlobalEigenVariations::initialize: dynamic cast failed\n";
        continue;
      }
      TH1* resultVariedUp   = (TH1*)result->Clone(eigenvarup);   resultVariedUp->SetDirectory(0); // copy flavour result, want to set bin contents according to the combined eigenvartion flavour block
      TH1* resultVariedDown = (TH1*)result->Clone(eigenvardown); resultVariedDown->SetDirectory(0);
      int up_to_bin = flav_bins[flavour];
      int current_bin = 1; // starting from 1 for filling histograms
      for(int flav_bin = bin_baseline+1 ; flav_bin < up_to_bin+1 ; flav_bin++){ // the +1's here are to deal with bin numbering yet again...
        Int_t bin = result->GetBin(1,current_bin); // increment along smaller (result) flavour variation TH1
        resultVariedUp->SetBinContent(bin, var.first->GetBinContent(flav_bin));     // Sets the result TH1 bin content to the eigenvariation bin content
        resultVariedDown->SetBinContent(bin, var.second->GetBinContent(flav_bin)); // In this way, we achieve flavour variations with proper binning
        current_bin+=1;
      }
      bin_baseline+=up_to_bin;
      m_flav_eigen[flavour].push_back(std::make_pair(resultVariedUp, resultVariedDown)); // <-------- add the flavour EV to the storage (indexed the same as in m_eigen!)
    }
  }
 
  for(const auto& f : m_flavours){
    if (m_verbose) std::cout << " " << f << " m_flav_eigen has " << m_flav_eigen[f].size() << std::endl;
  }
 
  m_initialized = true;
}

isExtrapolationVariation() [1/2]

bool CalibrationDataEigenVariations::isExtrapolationVariation ( unsigned int nameIndex ) const

inherited

flag whether the given index corresponds to an extrapolation variation

Definition at line 1103 of file CalibrationDataEigenVariations.cxx.

{
  // Verifies whether the given named variation index corresponds to the extrapolation
  // uncertainty.
 
  return (m_namedExtrapolation == int(nameIndex));
}

isExtrapolationVariation() [2/2]

bool CalibrationDataGlobalEigenVariations::isExtrapolationVariation	(	unsigned int	nameIndex,
		const std::string &	flavour ) const

Definition at line 2041 of file CalibrationDataEigenVariations.cxx.

{
  // Verifies whether the given named variation index corresponds to the extrapolation
  // uncertainty.
  int extrapIndex = m_flav_namedExtrapolation.find(flavour)->second;
  return (extrapIndex == int(nameIndex));
}

listNamedVariations() [1/2]

vector< string > CalibrationDataEigenVariations::listNamedVariations ( ) const

inherited

list the named variations

Definition at line 1001 of file CalibrationDataEigenVariations.cxx.

{
  // Provides the list of named variations
 
  vector<string> names;
  for (map<string, unsigned int>::const_iterator it = m_namedIndices.begin(); it != m_namedIndices.end(); ++it){
    names.push_back(it->first);
  }
  return names;
}

listNamedVariations() [2/2]

std::vector< std::string > CalibrationDataGlobalEigenVariations::listNamedVariations

(

const std::string &

flavour

)

const

Definition at line 1552 of file CalibrationDataEigenVariations.cxx.

{
  // Provides the list of named variations
 
  std::vector<std::string> names;
  for(const auto& namedar : m_flav_namedIndices.find(flavour)->second){
    names.push_back(namedar.first);
  }
  return names;
}

mergeVariations() [1/4]

void CalibrationDataEigenVariations::mergeVariations ( const IndexSet & set )

inherited

merge all variations in the given set

Definition at line 881 of file CalibrationDataEigenVariations.cxx.

{
  IndexSuperSet sset;
  sset.insert(set);
  mergeVariations(sset);
}

mergeVariations() [2/4]

void CalibrationDataEigenVariations::mergeVariations ( const IndexSuperSet & set )

inherited

merge all variations in any of the given sets

Definition at line 890 of file CalibrationDataEigenVariations.cxx.

{
  // check for overlap
  IndexSet checker;
  for (IndexSuperSet::iterator set_it = set.begin(); set_it != set.end(); ++set_it) {
    for (IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it){
      if (checker.count(*subset_it) == 0 && *subset_it <= m_eigen.size())
        checker.insert(*subset_it);
      else {
        std::cerr << "Error in CalibrationDataEigenVariations::mergeVariations: \
          IndexSets must not overlap and must lie between 1 and " << m_eigen.size() << ". Aborting!" << std::endl;
        return;
      }
    }
  }
 
  // retrieve the central calibration
  TH1 *result = static_cast<TH1*>(m_cnt->GetValue("result"));
  IndexSet toDelete;
  int nbins = result->GetNbinsX()+2;
  int ndim  = result->GetDimension();
  if (ndim > 1) nbins *= (result->GetNbinsY()+2);
  if (ndim > 2) nbins *= (result->GetNbinsZ()+2);
 
  // TH1 *var_up_final = static_cast<TH1*>(result->Clone()),
  //   *var_down_final = static_cast<TH1*>(result->Clone());
 
  // var_up_final->Reset();
  // var_down_final->Reset();
 
  // complex sum
  for (IndexSuperSet::iterator set_it = set.begin(); set_it != set.end(); ++set_it) {
    if (set_it->empty()) continue;
 
    double sum_H_up = 0.0, sum_H_down = 0.0;
    size_t lowest_index = *set_it->lower_bound(0);
    TH1 *total_var_up = static_cast<TH1*>(m_eigen[lowest_index].first->Clone()),
      *total_var_down = static_cast<TH1*>(m_eigen[lowest_index].second->Clone());
    total_var_up->SetDirectory(0);
    total_var_down->SetDirectory(0);
 
    total_var_up->Reset();
    total_var_down->Reset();
 
    // sum all other variations
    for (IndexSet::iterator subset_it = set_it->begin();
     subset_it != set_it->end(); ++subset_it) {
      size_t actual_index = *subset_it;
 
      if (actual_index != lowest_index) toDelete.insert(*subset_it);
 
      TH1 *partial_var_up = static_cast<TH1*>(m_eigen[actual_index].first->Clone()),
    *partial_var_down = static_cast<TH1*>(m_eigen[actual_index].second->Clone());
      partial_var_up->SetDirectory(0);
      partial_var_down->SetDirectory(0);
      
      partial_var_up->Add(result, -1.0); // <----- Is this correct? Should it be +1?
      partial_var_down->Add(result, -1.0);
      for (int i = 0; i < nbins; ++i) {
        partial_var_down->SetBinContent(i, -1.0*partial_var_down->GetBinContent(i));
      }
 
      for (int u = 0; u < nbins; ++u) {
        double sum_up = total_var_up->GetBinContent(u),
      sum_down = total_var_down->GetBinContent(u);
        for (int v = 0; v < nbins; ++v) {
          sum_up += partial_var_up->GetBinContent(u)*partial_var_up->GetBinContent(v);
          sum_H_up += partial_var_up->GetBinContent(u)*partial_var_up->GetBinContent(v);
          sum_down += partial_var_down->GetBinContent(u)*partial_var_down->GetBinContent(v);
          sum_H_down += partial_var_down->GetBinContent(u)*partial_var_down->GetBinContent(v);
        }
        total_var_up->SetBinContent(u, sum_up);
        total_var_down->SetBinContent(u, sum_down);
      }
      delete partial_var_up;
      delete partial_var_down;
    }
 
    // final part of complex summing
    for (int i = 0; i < nbins; ++i) {
      if (sum_H_up != 0.0)
        total_var_up->SetBinContent(i, total_var_up->GetBinContent(i)/sqrt(sum_H_up));
      else
        total_var_up->SetBinContent(i, 0.0);
      if (sum_H_down != 0.0)
        total_var_down->SetBinContent(i, -1.0*total_var_down->GetBinContent(i)/sqrt(sum_H_down));
      else
        total_var_down->SetBinContent(i, 0.0);
    }
 
    total_var_up->Add(result);
    total_var_down->Add(result);
 
    m_eigen[lowest_index].first = total_var_up;
    m_eigen[lowest_index].second = total_var_down;
  }
 
  removeVariations(toDelete);
}

mergeVariations() [3/4]

void CalibrationDataGlobalEigenVariations::mergeVariations	(	const IndexSet &	set,
		std::string &	flav )

Definition at line 2068 of file CalibrationDataEigenVariations.cxx.

{
  IndexSuperSet sset;
  sset.insert(set);
  mergeVariations(sset, flav);
}

mergeVariations() [4/4]

void CalibrationDataGlobalEigenVariations::mergeVariations	(	const IndexSuperSet &	set,
		std::string &	flav )

Definition at line 2078 of file CalibrationDataEigenVariations.cxx.

{
  // Merge the (flavour specific) eigenvariations, as stored in m_flav_eigen
 
  // check for overlap
  IndexSet checker;
  for (IndexSuperSet::iterator set_it = set.begin(); set_it != set.end(); ++set_it) {
    for (IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it){
      if (checker.count(*subset_it) == 0 && *subset_it <= m_flav_eigen[flav].size())
        checker.insert(*subset_it);
      else {
        std::cerr << "Error in CalibrationDataGlobalEigenVariations::mergeVariations: \
          IndexSets must not overlap and must lie between 1 and " << m_eigen.size() << ". Aborting!" << std::endl;
        return;
      }
    }
  }
 
  std::string flavour = flav;
  
  Analysis::CalibrationDataHistogramContainer* c = m_histcontainers[flavour];
 
  TH1* result = dynamic_cast<TH1*>(c->GetValue("result"));
  if (not result){
    std::cerr << "Error in CalibrationDataGlobalEigenVariations::mergeVariations: failed dynamic cast\n";
    return;
  }
  // retrieve the central calibration
 
  IndexSet toDelete;
  int nbins = result->GetNbinsX()+2;
  int ndim  = result->GetDimension();
  if (ndim > 1) nbins *= (result->GetNbinsY()+2);
  if (ndim > 2) nbins *= (result->GetNbinsZ()+2);
 
  // complex sum
  for (IndexSuperSet::iterator set_it = set.begin(); set_it != set.end(); ++set_it) {
    if (set_it->empty()) continue;
 
    double sum_H_up = 0.0, sum_H_down = 0.0;
    size_t lowest_index = *set_it->lower_bound(0);
    TH1 *total_var_up = static_cast<TH1*>(m_flav_eigen[flavour][lowest_index].first->Clone());
    TH1 *total_var_down = static_cast<TH1*>(m_flav_eigen[flavour][lowest_index].second->Clone());
    total_var_up->SetDirectory(0);
    total_var_down->SetDirectory(0);
 
    total_var_up->Reset();
    total_var_down->Reset();
 
    // sum all other variations
    for (IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it) {
      size_t actual_index = *subset_it;
 
      if (actual_index != lowest_index) toDelete.insert(*subset_it); // 
 
      TH1 *partial_var_up = static_cast<TH1*>(m_flav_eigen[flavour][actual_index].first->Clone()); 
      TH1 *partial_var_down = static_cast<TH1*>(m_flav_eigen[flavour][actual_index].second->Clone());
      partial_var_up->SetDirectory(0);
      partial_var_down->SetDirectory(0);
      
      partial_var_up->Add(result, -1.0);
      partial_var_down->Add(result, -1.0);
      for (int i = 0; i < nbins; ++i) {
        partial_var_down->SetBinContent(i, -1.0*partial_var_down->GetBinContent(i));
      }
 
      for (int u = 0; u < nbins; ++u) {
        double sum_up = total_var_up->GetBinContent(u);
        double sum_down = total_var_down->GetBinContent(u);
        for (int v = 0; v < nbins; ++v) {
          sum_up += partial_var_up->GetBinContent(u)*partial_var_up->GetBinContent(v);
          sum_H_up += partial_var_up->GetBinContent(u)*partial_var_up->GetBinContent(v);
          sum_down += partial_var_down->GetBinContent(u)*partial_var_down->GetBinContent(v);
          sum_H_down += partial_var_down->GetBinContent(u)*partial_var_down->GetBinContent(v);
        }
        total_var_up->SetBinContent(u, sum_up);
        total_var_down->SetBinContent(u, sum_down);
      }
      delete partial_var_up;
      delete partial_var_down;
    }
 
    // final part of complex summing
    for (int i = 0; i < nbins; ++i) {
      if (sum_H_up != 0.0)
        total_var_up->SetBinContent(i, total_var_up->GetBinContent(i)/sqrt(sum_H_up));
      else
        total_var_up->SetBinContent(i, 0.0);
      if (sum_H_down != 0.0)
        total_var_down->SetBinContent(i, -1.0*total_var_down->GetBinContent(i)/sqrt(sum_H_down));
      else
        total_var_down->SetBinContent(i, 0.0);
    }
 
    total_var_up->Add(result);
    total_var_down->Add(result);
 
    m_flav_eigen[flavour][lowest_index].first = total_var_up;
    m_flav_eigen[flavour][lowest_index].second = total_var_down;
  }
  
  removeVariations(toDelete, flavour);
}

mergeVariationsFrom() [1/2]

void CalibrationDataEigenVariations::mergeVariationsFrom ( const size_t & index )

inherited

merge all variations starting from the given index

Definition at line 867 of file CalibrationDataEigenVariations.cxx.

{
  // Merge all systematic variation starting from the given index.
  // The resulting merged variation replaces the first entry in the list
  // (i.e., the entry specified by the index).
  IndexSet simple_set;
 
  for (size_t it = index; it < m_eigen.size(); ++it)
    simple_set.insert(it);
  mergeVariations(simple_set);
}

mergeVariationsFrom() [2/2]

void CalibrationDataGlobalEigenVariations::mergeVariationsFrom	(	const size_t &	index,
		std::string &	flav )

Definition at line 2053 of file CalibrationDataEigenVariations.cxx.

{
  // Merge all systematic variation starting from the given index.
  // The resulting merged variation replaces the first entry in the list
  // (i.e., the entry specified by the index).
 
  IndexSet simple_set;
 
  for (size_t it = index; it < m_flav_eigen[flav].size(); ++it)
    simple_set.insert(it);
  mergeVariations(simple_set, flav);
}

removeVariations() [1/4]

void CalibrationDataEigenVariations::removeVariations ( const IndexSet & set )

inherited

remove all variations in the given set

Definition at line 839 of file CalibrationDataEigenVariations.cxx.

{
  if (set.size() == 0) return;
 
  std::vector<std::pair<TH1*, TH1*> > new_eigen;
  for (size_t index = 0; index < m_eigen.size(); ++index){
      if (set.count(index) == 0) new_eigen.push_back(m_eigen[index]);
      else { delete m_eigen[index].first; delete m_eigen[index].second; }
    }
  m_eigen = std::move(new_eigen);
}

removeVariations() [2/4]

void CalibrationDataEigenVariations::removeVariations ( const IndexSuperSet & set )

inherited

remove all variations in any of the given sets

Definition at line 853 of file CalibrationDataEigenVariations.cxx.

{
  IndexSet simple_set;
 
  for (IndexSuperSet::iterator set_it = set.begin(); set_it != set.end(); ++set_it) {
    for (IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it)
      simple_set.insert(*subset_it);
  }
 
  removeVariations(simple_set);
}

removeVariations() [3/4]

void CalibrationDataGlobalEigenVariations::removeVariations	(	const IndexSet &	set,
		std::string &	flav )

Definition at line 2187 of file CalibrationDataEigenVariations.cxx.

{
  if (set.size() == 0) return;
  std::vector<std::pair<TH1*, TH1*> > new_eigen;
  std::vector<std::pair<TH1*, TH1*>> eigen = m_flav_eigen[flav];
  for (size_t index = 0; index < eigen.size(); ++index){
      if (set.count(index) == 0) {
        new_eigen.push_back(eigen[index]);
      } else { 
        delete eigen[index].first; delete eigen[index].second; 
      }
    }
 
  m_flav_eigen[flav] = std::move(new_eigen);
 
}

removeVariations() [4/4]

void CalibrationDataGlobalEigenVariations::removeVariations	(	const IndexSuperSet &	set,
		std::string &	flav )

Definition at line 2206 of file CalibrationDataEigenVariations.cxx.

{
  IndexSet simple_set;
 
  for (IndexSuperSet::iterator set_it = set.begin(); set_it != set.end(); ++set_it) {
    for (IndexSet::iterator subset_it = set_it->begin(); subset_it != set_it->end(); ++subset_it)
      simple_set.insert(*subset_it);
  }
 
  removeVariations(simple_set, flav);
}

setVerbose()

void CalibrationDataEigenVariations::setVerbose ( bool verbose )

inherited

Definition at line 1255 of file CalibrationDataEigenVariations.cxx.

                                                            {
  m_verbose = verbose;
}

Member Data Documentation

m_all_shared_systematics

std::set<std::string> Analysis::CalibrationDataGlobalEigenVariations::m_all_shared_systematics

private

Definition at line 181 of file CalibrationDataEigenVariations.h.

m_blockmatrixsize

int Analysis::CalibrationDataGlobalEigenVariations::m_blockmatrixsize

private

Definition at line 178 of file CalibrationDataEigenVariations.h.

m_capturedvariance

double Analysis::CalibrationDataEigenVariations::m_capturedvariance

protectedinherited

Definition at line 137 of file CalibrationDataEigenVariations.h.

m_cdipath

std::string Analysis::CalibrationDataEigenVariations::m_cdipath

protectedinherited

@ data members needed for eigenvector method

the map stores the int which is needed to access the other vector<> objects

Definition at line 129 of file CalibrationDataEigenVariations.h.

m_cnt

CalibrationDataHistogramContainer* Analysis::CalibrationDataEigenVariations::m_cnt

privateinherited

container object containing the basic information

Definition at line 100 of file CalibrationDataEigenVariations.h.

m_eigen

std::vector<std::pair<TH1*, TH1*> > Analysis::CalibrationDataEigenVariations::m_eigen

protectedinherited

eigenvector variations

Definition at line 116 of file CalibrationDataEigenVariations.h.

m_flav_eigen

std::map<std::string, std::vector<std::pair<TH1*, TH1*> > > Analysis::CalibrationDataGlobalEigenVariations::m_flav_eigen

private

Definition at line 189 of file CalibrationDataEigenVariations.h.

m_flav_named

std::map<std::string, std::vector<std::pair<TH1*,TH1*> > > Analysis::CalibrationDataGlobalEigenVariations::m_flav_named

private

Definition at line 187 of file CalibrationDataEigenVariations.h.

m_flav_namedExtrapolation

std::map<std::string, int> Analysis::CalibrationDataGlobalEigenVariations::m_flav_namedExtrapolation

private

Definition at line 188 of file CalibrationDataEigenVariations.h.

m_flav_namedIndices

std::map<std::string, std::map<std::string, unsigned int> > Analysis::CalibrationDataGlobalEigenVariations::m_flav_namedIndices

private

Definition at line 186 of file CalibrationDataEigenVariations.h.

m_flavour_combinations

std::map<std::string, std::vector<int> > Analysis::CalibrationDataGlobalEigenVariations::m_flavour_combinations

private

Definition at line 183 of file CalibrationDataEigenVariations.h.

m_flavours

std::vector<std::string> Analysis::CalibrationDataGlobalEigenVariations::m_flavours

private

Definition at line 190 of file CalibrationDataEigenVariations.h.

m_histcontainers

std::map<std::string, Analysis::CalibrationDataHistogramContainer*> Analysis::CalibrationDataGlobalEigenVariations::m_histcontainers

private

Definition at line 180 of file CalibrationDataEigenVariations.h.

m_initialized

bool Analysis::CalibrationDataEigenVariations::m_initialized

protectedinherited

flag whether the initialization has been carried out

Definition at line 105 of file CalibrationDataEigenVariations.h.

m_jetauthor

std::string Analysis::CalibrationDataEigenVariations::m_jetauthor

protectedinherited

Definition at line 132 of file CalibrationDataEigenVariations.h.

m_named

std::vector<std::pair<TH1*, TH1*> > Analysis::CalibrationDataEigenVariations::m_named

protectedinherited

Definition at line 110 of file CalibrationDataEigenVariations.h.

m_namedExtrapolation

int Analysis::CalibrationDataEigenVariations::m_namedExtrapolation

protectedinherited

named variation index for the special case of extrapolation uncertainties

Definition at line 113 of file CalibrationDataEigenVariations.h.

m_namedIndices

std::map<std::string, unsigned int> Analysis::CalibrationDataEigenVariations::m_namedIndices

protectedinherited

named variations

Definition at line 109 of file CalibrationDataEigenVariations.h.

m_only_shared_systematics

std::set<std::string> Analysis::CalibrationDataGlobalEigenVariations::m_only_shared_systematics

private

Definition at line 182 of file CalibrationDataEigenVariations.h.

m_statVariations

bool Analysis::CalibrationDataEigenVariations::m_statVariations

protectedinherited

indicate whether statistical uncertainties are stored as variations

Definition at line 119 of file CalibrationDataEigenVariations.h.

m_taggername

std::string Analysis::CalibrationDataEigenVariations::m_taggername

protectedinherited

Definition at line 130 of file CalibrationDataEigenVariations.h.

m_totalvariance

double Analysis::CalibrationDataEigenVariations::m_totalvariance

protectedinherited

Definition at line 136 of file CalibrationDataEigenVariations.h.

m_validate

bool Analysis::CalibrationDataEigenVariations::m_validate

protectedinherited

Definition at line 106 of file CalibrationDataEigenVariations.h.

m_verbose

bool Analysis::CalibrationDataEigenVariations::m_verbose

protectedinherited

Definition at line 139 of file CalibrationDataEigenVariations.h.

m_wp

std::string Analysis::CalibrationDataEigenVariations::m_wp

protectedinherited

Definition at line 131 of file CalibrationDataEigenVariations.h.

---

The documentation for this class was generated from the following files:

• CalibrationDataEigenVariations.h
• CalibrationDataEigenVariations.cxx