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