ErrorPropagator.h

Go to the documentation of this file.

#ifndef AMPGEN_ERRORPROPAGATOR_H
#define AMPGEN_ERRORPROPAGATOR_H
 
#include <cmath>
#include <stddef.h>
#include <string>
#include <vector>
#include <map>
#include <functional>
#include <array>
#include <utility>
 
#include <TMatrixD.h>
#include <TRandom3.h>
 
#include "AmpGen/MinuitParameter.h"
#include "AmpGen/MsgService.h"
 
class TRandom3;
 
namespace AmpGen {
  class MinuitParameterSet;
  class LinearErrorPropagator {
  public:
    explicit LinearErrorPropagator(const std::vector<MinuitParameter *> &params);
    explicit LinearErrorPropagator(const MinuitParameterSet &params);
    LinearErrorPropagator(const TMatrixD &reducedCovarianceMatrix, const std::vector<MinuitParameter *> &params);

    template <class FCN> double operator()(FCN fcn) const { return getError(std::function<double(void)>(fcn)); }

    double getError(const std::function<double(void)> &fcn) const;

    std::vector<double> getVectorError(const std::function<std::vector<double>(void)> &fcn, size_t RANK) const;

    std::vector<double> combinationWeights(const std::vector<std::function<double(void)>> &functions);
    TMatrixD covarianceMatrix(const std::vector<std::function<double(void)>> &functions);

    TMatrixD correlationMatrix(const std::vector<std::function<double(void)>> &functions);

    std::pair<double, double> combinationCovWeighted(const std::vector<std::function<double(void)>> &functions);
 
    void add(const LinearErrorPropagator &p2);
    void reset();
    size_t size() const;
    const TMatrixD &cov() const;
    const std::map<std::string, size_t> posMap() const;
    const std::vector<MinuitParameter *> &params() const;
 
  private:
    TMatrixD m_cov;
    std::vector<MinuitParameter *> m_parameters;
 
    template <class FCN> double derivative(FCN fcn, const size_t &i) const {
      double startingValue = m_parameters[i]->mean();
      m_parameters[i]->setCurrentFitVal(startingValue + std::sqrt(m_cov(i, i)));
      double plus_variation = fcn();
      m_parameters[i]->setCurrentFitVal(startingValue - std::sqrt(m_cov(i, i)));
      double minus_variation = fcn();
      m_parameters[i]->setCurrentFitVal(startingValue);
      return (plus_variation - minus_variation) / (2 * std::sqrt(m_cov(i, i)));
    }
  };
 
  class GaussErrorPropagator {
  public:
    GaussErrorPropagator(const TMatrixD &reducedCovariance, const std::vector<MinuitParameter *> &params, TRandom3 *rnd);
    void perturb();
    void reset();
    void transpose();
 
  private:
    std::vector<MinuitParameter *> m_parameters;
    std::vector<double> m_startingValues;
    TRandom3 *m_rand;
    TMatrixD m_decomposedCholesky;
  };
 
  class NonlinearErrorPropagator {
  public:
    NonlinearErrorPropagator(const std::function<double(void)> &, const TMatrixD &, const std::vector<MinuitParameter *> &);
    TMatrixD correlationMatrix(const std::vector<std::function<double(void)>> &functions, TRandom3 *rnd, const unsigned &nSamples = 50000);
 
  private:
    std::function<double(void)> m_fcn;
    TMatrixD m_cov;
    std::vector<MinuitParameter *> m_parameters;
  };
 
} // namespace AmpGen
#endif