UnfoldTikhonov.cxx
Go to the documentation of this file.
2 
4 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/Var.h"
7 
8 #include "Minuit2/MnMigrad.h"
9 #include "Minuit2/FunctionMinimum.h"
10 #include "Minuit2/MnUserParameters.h"
11 
12 #include "TError.h"
13 #include "TH2.h"
14 
15 #include <iostream>
16 
17 namespace ana
18 {
19  //----------------------------------------------------------------------
21  double reg)
22  : fNorm(-1), fReg(reg)
23  {
24  DontAddDirectory guard;
25 
26  std::unique_ptr<TH2D> hrt(rt.ToTH2(1e20)); // normalization doesn't matter
27 
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);
33  }
34  }
35 
36  fTrueBinning = hrt->ProjectionY(UniqueName().c_str());
37  fTrueBinning->Reset();
38 
39  if(fReg < 0){
41  std::cout << "Automatic regularization strength = " << fReg << std::endl;
42  }
43  }
44 
45  //----------------------------------------------------------------------
47  {
48  delete fTrueBinning;
49  }
50 
51  //----------------------------------------------------------------------
53  {
54  // Don't scale anything, ie truth prediction is MC truth
55  const std::vector<double> pars(NTrueBins(), 1);
56  std::vector<double> reco, truth;
57  RecoTruthPrediction(pars, reco, truth);
58 
59  // What would be the penalty at that prediction
60  const double pen = PenaltyTerm(truth);
61 
62  // Statistical fluctuations imply chisq = NDOF
63  const int chisq = NRecoBins();
64 
65  // Set coefficient so penalty term and chisq match
66  return chisq/pen;
67  }
68 
69  //----------------------------------------------------------------------
71  {
72  // Parameters are weights to apply to the initial true distribution
73  ROOT::Minuit2::MnUserParameters mnPars;
74  for(int i = 0; i < NTrueBins(); ++i){
75  mnPars.Add(TString::Format("w%d", i).Data(), 1, .5);
76  }
77 
78  std::unique_ptr<TH1D> hreco(reco.ToTH1(reco.POT()));
79  assert(hreco->GetNbinsX() == NRecoBins());
80  fReco.clear();
81  fReco.resize(NRecoBins());
82  fNorm = 0;
83  for(int i = 0; i < NRecoBins(); ++i){
84  fReco[i] = hreco->GetBinContent(i+1);
85  fNorm += fReco[i];
86  }
87 
88  ROOT::Minuit2::MnMigrad mnMigrad(*this, mnPars);
89  const int olderr = gErrorIgnoreLevel;
90  gErrorIgnoreLevel = 1001; // Ignore warnings
91  ROOT::Minuit2::FunctionMinimum minpt = mnMigrad();
92  gErrorIgnoreLevel = olderr;
93 
94  const std::vector<double> minvec = minpt.UserParameters().Params();
95 
96  std::vector<double> pred_reco, pred_truth;
97  RecoTruthPrediction(minvec, pred_reco, pred_truth);
98 
99  fReco.clear();
100  fNorm = -1;
101 
102  Eigen::ArrayXd ret(pred_truth.size()+2);
103  ret[0] = 0;
104  ret[pred_truth.size()+1] = 0;
105  for(unsigned int i = 0; i < pred_truth.size(); ++i) ret[i+1] = pred_truth[i];
106 
107  return Spectrum(ret,
108  HistAxis(reco.GetLabels(), reco.GetBinnings()),
109  reco.POT(), reco.Livetime());
110  }
111 
112  //----------------------------------------------------------------------
113  void UnfoldTikhonov::RecoTruthPrediction(const std::vector<double>& pars,
114  std::vector<double>& reco,
115  std::vector<double>& truth) const
116  {
117  const int X = NRecoBins();
118  const int Y = NTrueBins();
119 
120  reco.resize(X, 0);
121  truth.resize(Y, 0);
122 
123  double ws[Y];
124  for(int iy = 0; iy < Y; ++iy) ws[iy] = std::max(pars[iy], 0.);
125 
126  double tot = 0;
127 
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;
132  truth[iy] += val;
133  reco[ix] += val;
134  tot += val;
135  }
136  }
137 
138  // Ensure result normalization
139  if(fNorm > 0){
140  for(int ix = 0; ix < X; ++ix) reco [ix] *= fNorm/tot;
141  for(int iy = 0; iy < Y; ++iy) truth[iy] *= fNorm/tot;
142  }
143  }
144 
145  //----------------------------------------------------------------------
146  double UnfoldTikhonov::PenaltyTerm(const std::vector<double>& pred) const
147  {
148  double pen = 0;
149 
150  for(unsigned int i = 1; i+1 < pred.size(); ++i){
151  // Second derivative
152  const double diff2 = pred[i+1] -2*pred[i] + pred[i-1];
153 
154  // Normalize to absolute function value
155  const double val = pred[i];
156 
157  // Sum of squares
158  if(val) pen += diff2*diff2 / (val*val);
159  }
160 
161  return pen;
162  }
163 
164  //----------------------------------------------------------------------
165  double UnfoldTikhonov::ChiSq(const std::vector<double>& pred) const
166  {
167  double chi = 0;
168 
169  // Discourage negative values
170  for(unsigned int i = 0; i < pred.size(); ++i)
171  if(pred[i] < 0) chi += pred[i]*pred[i];
172 
173  for(int ix = 0; ix < NRecoBins(); ++ix){
174  const double dat = fReco[ix];
175 
176  double tot = pred[ix];
177 
178  chi += 2*(tot-dat);
179  if(dat) chi += 2*dat*log(dat/tot);
180  }
181 
182  return chi;
183  }
184 
185  //----------------------------------------------------------------------
186  double UnfoldTikhonov::operator()(const std::vector<double>& pars) const
187  {
188  // static int N = 0;
189  // if(N++%1000 == 0) std::cout << "." << std::flush;
190 
191  assert(int(pars.size()) == NTrueBins());
192 
193  std::vector<double> reco_pred, truth_pred;
194  RecoTruthPrediction(pars, reco_pred, truth_pred);
195 
196  const double chi_stat = ChiSq(reco_pred);
197  const double pen = fReg * PenaltyTerm(truth_pred);
198 
199  return chi_stat + pen;
200  }
201 }
T max(const caf::Proxy< T > &a, T b)
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:256
int NTrueBins() const
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
double ChiSq(const std::vector< double > &pred) const
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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.
Definition: Spectrum.cxx:148
UnfoldTikhonov(const ReweightableSpectrum &rt, double reg=-1)
Spectrum with the value of a second variable, allowing for reweighting
void RecoTruthPrediction(const std::vector< double > &pars, std::vector< double > &reco, std::vector< double > &truth) const
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Float_t Y
Definition: plot.C:38
double fReg
Regularization strength.
std::vector< double > fReco
std::vector< float > Spectrum
Definition: Constants.h:570
double POT() const
Definition: Spectrum.h:219
virtual Spectrum Truth(const Spectrum &reco) override
OStream cout
Definition: OStream.cxx:6
int NRecoBins() const
virtual double PenaltyTerm(const std::vector< double > &pred) const
std::string pars("Th23Dmsq32")
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)
Definition: Style.cxx:154
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:255
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t X
Definition: plot.C:38
Float_t w
Definition: plot.C:20
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:222
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
std::vector< std::vector< double > > fRT
TH2D * ToTH2(double pot) const
virtual double operator()(const std::vector< double > &pars) const override