10 #include "TDirectory.h" 11 #include "TObjString.h" 19 #include <Eigen/Dense> 30 using std::runtime_error;
31 using std::setprecision;
36 using std::chrono::high_resolution_clock;
37 using std::chrono::duration_cast;
38 using std::chrono::milliseconds;
39 using std::chrono::seconds;
40 using std::chrono::minutes;
44 using Eigen::MatrixXd;
45 using Eigen::VectorXd;
50 using covmx::CovarianceMatrix;
54 CovarianceMatrix* covmx,
double epsilon,
double lambdazero,
double nu)
55 : fSamples(samples), fCovMx(covmx), fNBins(0), fNBinsFull(0), fLambdaZero(lambdazero),
56 fNu(nu), fBetaHist(nullptr), fMuHist(nullptr),
57 fBetaMuHist(nullptr), fPredShiftedHist(nullptr), fPredUnshiftedHist(nullptr),
58 fDataHist(nullptr), fResetShifts(true), fVerbose(false)
61 fNBins +=
s.GetBinning().NBins();
82 for (
size_t iS = 0; iS < fSamples.size(); ++iS) {
83 double pot = fSamples[iS].GetData()->POT();
84 double livetime = fSamples[iS].GetData()->Livetime();
85 TH1D* hData = fSamples[iS].GetData()->ToTH1(pot);
86 TH1D* hCos = fSamples[iS].HasCosmic() ?
87 fSamples[iS].GetCosmic()->ToTH1(livetime,
kLivetime) :
nullptr;
89 fData(i) = hData->GetBinContent(iBin);
90 fCosmic(i) = hCos ? hCos->GetBinContent(iBin) : 0;
117 for (
size_t iS = 0; iS <
fSamples.size(); ++iS) {
122 shift, get<0>(comp), get<1>(comp), get<2>(comp));
123 TH1D* hMu = spec.
ToTH1(pot);
124 for (
size_t iBin = 1; iBin <= (size_t)hMu->GetNbinsX(); ++iBin) {
125 fMu(i) = hMu->GetBinContent(iBin);
134 if (calc &&
fMu.array().isNaN().any()) {
135 cout <<
"ERROR! NAN VALUES IN PREDICTION!" <<
endl;
136 cout <<
"---------- OSCILLATION PARAMETERS ----------" << endl
137 <<
"Dm21: " << calc->
GetDm(2) << endl
138 <<
"Dm31: " << calc->
GetDm(3) << endl
139 <<
"Dm41: " << calc->
GetDm(4) << endl
140 <<
"Theta 13: " << calc->
GetAngle(1, 3) << endl
141 <<
"Theta 23: " << calc->
GetAngle(2, 3) << endl
142 <<
"Theta 14: " << calc->
GetAngle(1, 4) << endl
143 <<
"Theta 24: " << calc->
GetAngle(2, 4) << endl
144 <<
"Theta 34: " << calc->
GetAngle(3, 4) << endl
145 <<
"Delta 13: " << calc->
GetDelta(1, 3) << endl
146 <<
"Delta 24: " << calc->
GetDelta(2, 4) << endl
147 <<
"--------------------------------------------" <<
endl;
163 "Number of bins in predictions does not match covariance matrix!");
182 unsigned int fullOffset(0), scaledOffset(0);
184 for (
size_t iS = 0; iS <
fSamples.size(); ++iS) {
185 for (
size_t iC = 0; iC <
fComps[iS].size(); ++iC) {
187 unsigned int iScaled = scaledOffset + iB;
188 unsigned int iFull = fullOffset + (iC*fNBinsPerSample[iS]) + iB;
195 unsigned int iScaled = scaledOffset + iB;
198 scaledOffset += fNBinsPerSample[iS];
199 fullOffset +=
fComps[iS].size() * fNBinsPerSample[iS];
221 unsigned int fullOffset(0), scaledOffset(0);
222 for (
size_t iS = 0; iS <
fSamples.size(); ++iS) {
224 unsigned int iScaled = scaledOffset+iB;
226 for (
size_t iC = 0; iC <
fComps[iS].size(); ++iC) {
227 unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
230 double val =
fData(iScaled) / (sum + eps);
231 for (
size_t iC = 0; iC <
fComps[iS].size(); ++iC) {
232 unsigned int iFull = fullOffset+(iC*fNBinsPerSample[iS])+iB;
236 scaledOffset += fNBinsPerSample[iS];
237 fullOffset +=
fComps[iS].size() * fNBinsPerSample[iS];
245 bool maskChange =
false;
247 size_t iFullOffset(0.), iScaledOffset(0.);
248 for (
size_t iS = 0; iS <
fSamples.size(); ++iS) {
249 for (
size_t iC = 0; iC <
fComps[iS].size(); ++iC) {
252 unsigned int iScaled = iScaledOffset+iB;
253 unsigned int iFull = iFullOffset+(iC*fNBinsPerSample[iS])+iB;
256 for (
size_t jC = 0; jC <
fComps[iS].size(); ++jC) {
258 unsigned int jFull = iFullOffset+(jC*fNBinsPerSample[iS])+iB;
262 if (y < 1
e-40) y = 1
e-40;
265 for (
size_t jFull = 0; jFull <
fNBinsFull; ++jFull) {
266 double beta = (jFull == iFull) ? 0 :
fBeta[jFull];
267 z += (beta-1)*
fMxInv(iFull, jFull);
271 double gradZero = 2 * (
fMu(iFull)*(1.0-(
fData(iScaled)/
y)) +
z);
276 if (gradZero <= 0 && !
fBetaMask[iFull]) {
294 size_t iFullOffset(0.), iScaledOffset(0.);
295 for (
size_t iS = 0; iS <
fSamples.size(); ++iS) {
296 for (
size_t iC = 0; iC <
fComps[iS].size(); ++iC) {
299 unsigned int iScaled = iScaledOffset+iB;
300 unsigned int iFull = iFullOffset+(iC*fNBinsPerSample[iS])+iB;
302 for (
size_t jC = 0; jC <
fComps[iS].size(); ++jC) {
303 unsigned int jFull = iFullOffset+(jC*fNBinsPerSample[iS])+iB;
306 if (y < 1
e-40) y = 1
e-40;
308 for (
size_t jFull = 0; jFull <
fNBinsFull; ++jFull) {
316 size_t jFullOffset(0.), jScaledOffset(0.);
317 for (
size_t jS = 0; jS <
fSamples.size(); ++jS) {
318 for (
size_t jC = 0; jC <
fComps[jS].size(); ++jC) {
319 for (
size_t jB = 0; jB < fNBinsPerSample[jS]; ++jB) {
321 unsigned int jScaled = jScaledOffset + jB;
322 unsigned int jFull = jFullOffset+(jC*fNBinsPerSample[jS])+jB;
323 if (iScaled != jScaled) {
331 jScaledOffset += fNBinsPerSample[jS];
332 jFullOffset +=
fComps[jS].size() * fNBinsPerSample[jS];
396 const int maxIters = 5e3;
406 bool keepChangingMask =
true;
407 vector<vector<bool>> maskHistory;
408 vector<double> chisqHistory;
412 double minChi =
prev;
414 cout <<
"Before minimization, chisq is " << prev
422 cout <<
"No convergence after " << maxIters <<
" iterations! Quitting out of the infinite loop." <<
endl;
429 if (deltaBeta.array().isNaN().any()) {
431 for (
int i = 0;
i < deltaBeta.size(); ++
i) {
432 if (
isnan(deltaBeta[
i])) ++nanCounter;
437 cout <<
"LDLT solution failed! There are " << nanCounter <<
"nan delta betas. Resetting all betas.";
452 fBeta(i) += deltaBeta(counter);
456 }
else if (
fBeta(i) > 1e10) {
470 if (
isnan(chisq) || chisq > 1e20) {
472 cout <<
"ChiSq is NaN! Resetting minimisation." <<
endl;
473 else if (
fVerbose && chisq > 1e20) {
474 cout <<
"ChiSq is anomalously large! Resetting minimisation." <<
endl;
475 vector<double>
changes(deltaBeta.size());
476 VectorXd::Map(&
changes[0], deltaBeta.size()) = deltaBeta;
477 std::sort(
changes.begin(),
changes.end(), std::greater<double>());
478 cout <<
" five greatest deltas:";
484 fLambda *= (1 + ((double)nFailures/10.));
491 if (chisq < minChi) minChi = chisq;
502 assert(
false &&
"Negative systematic penalty!");
506 double change =
fabs(chisq-prev);
509 if (change/chisq < 1e-3 || change < 1e-10) {
511 bool maskChange = keepChangingMask ?
MaskBetas() :
false;
517 chisqHistory.push_back(chisq);
518 for (
size_t i = 0;
i < maskHistory.size()-1; ++
i) {
520 if (
fVerbose)
cout <<
"Found a repeated mask! Setting the mask that provided the best chisq" <<
endl;
523 for (
size_t j = 0;
j < chisqHistory.size(); ++
j) {
524 if (chisqHistory[
j] < chi) {
525 chi = chisqHistory[
j];
531 keepChangingMask =
false;
538 cout <<
"Minimisation round finished with chisq " << chisq <<
", bin mask:";
546 if (!maskChange && (change/chisq < 1e-10 || change < 1e-10)) {
550 <<
" syst) after " <<
fIteration <<
" iterations and ";
554 double millisecs = duration_cast<milliseconds>(
end-
start).
count();
555 cout << millisecs <<
" ms." <<
endl;
556 }
else if (secs < 60)
cout << secs <<
" seconds." <<
endl;
559 cout << mins <<
" minutes." <<
endl;
621 fDataHist->SetBinContent(
i+1,
fData(
i));
double fLambda
Levenberg-Marquardt nu.
Eigen::MatrixXd fHessReduced
double GetDelta(int i, int j) const
void InitialiseBetas() const
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
std::vector< std::vector< covmx::Component > > fComps
double LikelihoodCovMxNewton() const
LikelihoodCovMxExperiment(std::vector< covmx::Sample > samples, covmx::CovarianceMatrix *covmx, double epsilon=1e-5, double lambdazero=0, double nu=10)
void IncreaseLambda() const
fvar< T > fabs(const fvar< T > &x)
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
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.
double fNu
Levenberg-Marquardt starting lambda.
Eigen::MatrixXd GetFullCovMx()
covmx::CovarianceMatrix * fCovMx
Simple record of shifts applied to systematic parameters.
bool equal(T *first, T *second)
void DecreaseLambda() const
Adapt the PMNS_Sterile calculator to standard interface.
TH1D * fBetaHist
Levenberg-Marquardt lambda.
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
std::tuple< Flavors::Flavors_t, Current::Current_t, Sign::Sign_t > Component
Representation of a spectrum in any variable, with associated POT.
int isnan(const stan::math::var &a)
double GetChiSq(Eigen::ArrayXd e) const
TH1D * fPredUnshiftedHist
void GetReducedGradAndHess() const
Oscillation probability calculators.
double fLambdaZero
Whether to use LMA.
const std::vector< double > & Edges() const
void SaveHists(bool opt=true)
std::vector< unsigned int > fNBinsPerSample
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
void GetGradAndHess() const
std::vector< covmx::Sample > fSamples
assert(nhit_max >=nhit_nbins)
double GetDm(int i) const
std::vector< bool > fBetaMask
Eigen::ArrayXd GetExpectedSpectrum() const
Prevent histograms being added to the current directory.
std::string UniqueName()
Return a different string each time, for creating histograms.
double GetAngle(int i, int j) const
Eigen::VectorXd fGradReduced