Teddy.h

Go to the documentation of this file.

#pragma once
 
#include <functional>
#include <algorithm>
#include <memory>
 
#include <Eigen/Dense>
 
#include <TRandom3.h>
#include <TH2.h>
 
class TUnfoldBinning;
 
namespace DAS {
 
Eigen::VectorXd identity (const Eigen::VectorXd& x);
 
class Teddy {
 
    static Eigen::VectorXd H2V (const std::unique_ptr<TH1D>&); 
    static Eigen::MatrixXd H2M (const std::unique_ptr<TH2D>&); 
 
    static int GuessOutDim (int, std::function<Eigen::VectorXd(const Eigen::VectorXd&)>);
 
    static void SanityCheck (const Eigen::VectorXd&, Eigen::VectorXd&, double = 0);
 
    const Eigen::VectorXd vec; 
    const Eigen::MatrixXd cov; 
 
    const int nOldBins, 
              nNewBins; 
    long long N; 
 
    Eigen::VectorXd eigVals; 
    Eigen::MatrixXd rotation; 
 
    Eigen::VectorXd sigma; 
 
    Eigen::VectorXd sum; 
    Eigen::MatrixXd sum2; 
 
    std::function<Eigen::VectorXd(const Eigen::VectorXd &)> func; 
 
    TRandom3 random; 
 
public:
 
    static bool verbose;
 
    Teddy
        (const Eigen::VectorXd&, 
         const Eigen::MatrixXd&, 
         std::function<Eigen::VectorXd(const Eigen::VectorXd &)> = identity, 
         int = 42 
         );
 
    Teddy
        (const std::unique_ptr<TH1D>&, 
         const std::unique_ptr<TH2D>&, 
         std::function<Eigen::VectorXd(const Eigen::VectorXd &)> = identity, 
         int = 42 
         );
 
    void play (); 
    std::pair<Eigen::VectorXd, Eigen::MatrixXd> buy () const; 
 
    template<typename... Args>
        std::pair<std::unique_ptr<TH1D>,std::unique_ptr<TH2D>> buy 
            (TString name, Args... args) const
    {
        using namespace std;
        using namespace Eigen;
 
        TAxis * axis = new TAxis(args...);
        int ndiv = axis->GetNbins();
        double * edges = new double[ndiv];
        axis->GetLowEdge(edges);
        edges[ndiv-1] = axis->GetBinUpEdge(ndiv);
        auto h   = make_unique<TH1D>(name         , name, ndiv-1, edges);
        auto cov = make_unique<TH2D>(name + "_cov", name, ndiv-1, edges, ndiv-1, edges);
 
        VectorXd x(nNewBins);
        MatrixXd err(nNewBins, nNewBins);
        tie(x,err) = buy();
    
        for (int i=0; i<nNewBins; ++i) {
            double content = x(i),
                   error = sqrt(err(i,i));
    
            h->SetBinContent(i+1, content);
            h->SetBinError  (i+1, error  );
    
            for (int j=0; j<nNewBins; ++j)
                cov->SetBinContent(i+1,j+1, err(i,j));
        }
 
        return make_pair<unique_ptr<TH1D>,unique_ptr<TH2D>>(move(h),move(cov));
    }
};
bool Teddy::verbose = true;
 
} // end of DAS namespace