UnfoldTikhonov.cxx

Go to the documentation of this file.

#include "CAFAna/Unfold/UnfoldTikhonov.h"

#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/Var.h"

#include "Minuit2/MnMigrad.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameters.h"

#include "TError.h"
#include "TH2.h"

#include <iostream>

namespace ana
{
  //----------------------------------------------------------------------
  UnfoldTikhonov::UnfoldTikhonov(const ReweightableSpectrum& rt,
                                 double reg)
    : fNorm(-1), fReg(reg)
  {
    DontAddDirectory guard;

    std::unique_ptr<TH2D> hrt(rt.ToTH2(1e20)); // normalization doesn't matter

    fRT.resize(hrt->GetNbinsX());
    for(int ix = 0; ix < hrt->GetNbinsX(); ++ix){
      fRT[ix].resize(hrt->GetNbinsY());
      for(int iy = 0; iy < hrt->GetNbinsY(); ++iy){
        fRT[ix][iy] = hrt->GetBinContent(ix+1, iy+1);
      }
    }

    fTrueBinning = hrt->ProjectionY(UniqueName().c_str());
    fTrueBinning->Reset();

    if(fReg < 0){
      fReg = DeriveCrudeRegStrength();
      std::cout << "Automatic regularization strength = " << fReg << std::endl;
    }
  }

  //----------------------------------------------------------------------
  UnfoldTikhonov::~UnfoldTikhonov()
  {
    delete fTrueBinning;
  }

  //----------------------------------------------------------------------
  double UnfoldTikhonov::DeriveCrudeRegStrength() const
  {
    // Don't scale anything, ie truth prediction is MC truth
    const std::vector<double> pars(NTrueBins(), 1);
    std::vector<double> reco, truth;
    RecoTruthPrediction(pars, reco, truth);

    // What would be the penalty at that prediction
    const double pen = PenaltyTerm(truth);

    // Statistical fluctuations imply chisq = NDOF
    const int chisq = NRecoBins();

    // Set coefficient so penalty term and chisq match
    return chisq/pen;
  }

  //----------------------------------------------------------------------
  Spectrum UnfoldTikhonov::Truth(const Spectrum& reco)
  {
    // Parameters are weights to apply to the initial true distribution
    ROOT::Minuit2::MnUserParameters mnPars;
    for(int i = 0; i < NTrueBins(); ++i){
      mnPars.Add(TString::Format("w%d", i).Data(), 1, .5);
    }

    std::unique_ptr<TH1D> hreco(reco.ToTH1(reco.POT()));
    assert(hreco->GetNbinsX() == NRecoBins());
    fReco.clear();
    fReco.resize(NRecoBins());
    fNorm = 0;
    for(int i = 0; i < NRecoBins(); ++i){
      fReco[i] = hreco->GetBinContent(i+1);
      fNorm += fReco[i];
    }

    ROOT::Minuit2::MnMigrad mnMigrad(*this, mnPars);
    const int olderr = gErrorIgnoreLevel;
    gErrorIgnoreLevel = 1001; // Ignore warnings
    ROOT::Minuit2::FunctionMinimum minpt = mnMigrad();
    gErrorIgnoreLevel = olderr;

    const std::vector<double> minvec = minpt.UserParameters().Params();

    std::vector<double> pred_reco, pred_truth;
    RecoTruthPrediction(minvec, pred_reco, pred_truth);

    fReco.clear();
    fNorm = -1;

    Eigen::ArrayXd ret(pred_truth.size()+2);
    ret[0] = 0;
    ret[pred_truth.size()+1] = 0;
    for(unsigned int i = 0; i < pred_truth.size(); ++i) ret[i+1] = pred_truth[i];

    return Spectrum(ret,
                    HistAxis(reco.GetLabels(), reco.GetBinnings()),
                    reco.POT(), reco.Livetime());
  }

  //----------------------------------------------------------------------
  void UnfoldTikhonov::RecoTruthPrediction(const std::vector<double>& pars,
                                           std::vector<double>& reco,
                                           std::vector<double>& truth) const
  {
    const int X = NRecoBins();
    const int Y = NTrueBins();

    reco.resize(X, 0);
    truth.resize(Y, 0);

    double ws[Y];
    for(int iy = 0; iy < Y; ++iy) ws[iy] = std::max(pars[iy], 0.);

    double tot = 0;

    for(int iy = 0; iy < Y; ++iy){
      const double w = ws[iy];
      for(int ix = 0; ix < X; ++ix){
        const double val = fRT[ix][iy] * w;
        truth[iy] += val;
        reco[ix] += val;
        tot += val;
      }
    }

    // Ensure result normalization
    if(fNorm > 0){
      for(int ix = 0; ix < X; ++ix) reco [ix] *= fNorm/tot;
      for(int iy = 0; iy < Y; ++iy) truth[iy] *= fNorm/tot;
    }
  }

  //----------------------------------------------------------------------
  double UnfoldTikhonov::PenaltyTerm(const std::vector<double>& pred) const
  {
    double pen = 0;

    for(unsigned int i = 1; i+1 < pred.size(); ++i){
      // Second derivative
      const double diff2 = pred[i+1] -2*pred[i] + pred[i-1];

      // Normalize to absolute function value
      const double val = pred[i];

      // Sum of squares
      if(val) pen += diff2*diff2 / (val*val);
    }

    return pen;
  }

  //----------------------------------------------------------------------
  double UnfoldTikhonov::ChiSq(const std::vector<double>& pred) const
  {
    double chi = 0;

    // Discourage negative values
    for(unsigned int i = 0; i < pred.size(); ++i)
      if(pred[i] < 0) chi += pred[i]*pred[i];

    for(int ix = 0; ix < NRecoBins(); ++ix){
      const double dat = fReco[ix];

      double tot = pred[ix];

      chi += 2*(tot-dat);
      if(dat) chi += 2*dat*log(dat/tot);
    }

    return chi;
  }

  //----------------------------------------------------------------------
  double UnfoldTikhonov::operator()(const std::vector<double>& pars) const
  {
    //    static int N = 0;
    //    if(N++%1000 == 0) std::cout << "." << std::flush;

    assert(int(pars.size()) == NTrueBins());

    std::vector<double> reco_pred, truth_pred;
    RecoTruthPrediction(pars, reco_pred, truth_pred);

    const double chi_stat = ChiSq(reco_pred);
    const double pen = fReg * PenaltyTerm(truth_pred);

    return chi_stat + pen;
  }
}