6 #include "CAFAna/Core/Var.h" 8 #include "Minuit2/MnMigrad.h" 9 #include "Minuit2/FunctionMinimum.h" 10 #include "Minuit2/MnUserParameters.h" 22 : fNorm(-1), fReg(reg)
26 std::unique_ptr<TH2D> hrt(rt.
ToTH2(1e20));
28 fRT.resize(hrt->GetNbinsX());
29 for(
int ix = 0; ix < hrt->GetNbinsX(); ++ix){
30 fRT[ix].resize(hrt->GetNbinsY());
31 for(
int iy = 0; iy < hrt->GetNbinsY(); ++iy){
32 fRT[ix][iy] = hrt->GetBinContent(ix+1, iy+1);
56 std::vector<double> reco, truth;
73 ROOT::Minuit2::MnUserParameters mnPars;
78 std::unique_ptr<TH1D> hreco(reco.
ToTH1(reco.
POT()));
84 fReco[
i] = hreco->GetBinContent(
i+1);
88 ROOT::Minuit2::MnMigrad mnMigrad(*
this, mnPars);
91 ROOT::Minuit2::FunctionMinimum minpt = mnMigrad();
94 const std::vector<double> minvec = minpt.UserParameters().Params();
96 std::vector<double> pred_reco, pred_truth;
102 Eigen::ArrayXd
ret(pred_truth.size()+2);
104 ret[pred_truth.size()+1] = 0;
105 for(
unsigned int i = 0;
i < pred_truth.size(); ++
i) ret[
i+1] = pred_truth[
i];
114 std::vector<double>& reco,
115 std::vector<double>& truth)
const 124 for(
int iy = 0; iy <
Y; ++iy) ws[iy] =
std::max(pars[iy], 0.);
128 for(
int iy = 0; iy <
Y; ++iy){
129 const double w = ws[iy];
130 for(
int ix = 0; ix <
X; ++ix){
131 const double val =
fRT[ix][iy] *
w;
140 for(
int ix = 0; ix <
X; ++ix) reco [ix] *=
fNorm/tot;
141 for(
int iy = 0; iy <
Y; ++iy) truth[iy] *=
fNorm/tot;
150 for(
unsigned int i = 1;
i+1 < pred.size(); ++
i){
152 const double diff2 = pred[
i+1] -2*pred[
i] + pred[
i-1];
155 const double val = pred[
i];
158 if(val) pen += diff2*diff2 / (val*
val);
170 for(
unsigned int i = 0;
i < pred.size(); ++
i)
171 if(pred[
i] < 0) chi += pred[
i]*pred[
i];
176 double tot = pred[ix];
179 if(dat) chi += 2*dat*
log(dat/tot);
193 std::vector<double> reco_pred, truth_pred;
196 const double chi_stat =
ChiSq(reco_pred);
199 return chi_stat + pen;
T max(const caf::Proxy< T > &a, T b)
const std::vector< Binning > & GetBinnings() const
double ChiSq(const std::vector< double > &pred) const
Cuts and Vars for the 2020 FD DiF Study.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
UnfoldTikhonov(const ReweightableSpectrum &rt, double reg=-1)
Spectrum with the value of a second variable, allowing for reweighting
std::vector< double > Spectrum
void RecoTruthPrediction(const std::vector< double > &pars, std::vector< double > &reco, std::vector< double > &truth) const
virtual ~UnfoldTikhonov()
Representation of a spectrum in any variable, with associated POT.
double fReg
Regularization strength.
std::vector< double > fReco
virtual Spectrum Truth(const Spectrum &reco) override
virtual double PenaltyTerm(const std::vector< double > &pred) const
double fNorm
Sum of fReco.
assert(nhit_max >=nhit_nbins)
double DeriveCrudeRegStrength() const
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const std::vector< std::string > & GetLabels() const
Prevent histograms being added to the current directory.
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
_HistAxis< Var > HistAxis
std::string UniqueName()
Return a different string each time, for creating histograms.
std::vector< std::vector< double > > fRT
TH2D * ToTH2(double pot) const
virtual double operator()(const std::vector< double > &pars) const override