15 #include <Eigen/Dense> 23 using std::back_inserter;
30 using std::make_unique;
31 using std::ostringstream;
32 using std::runtime_error;
36 using std::unique_ptr;
39 using Eigen::MatrixXd;
41 using covmx::CovarianceMatrix;
50 transform(begin(samples),
end(samples), back_inserter(fSampleIDs),
51 [](
Sample s) ->
unsigned int {
return s.GetID(); });
52 if (!dir)
throw runtime_error(
"Directory passed is not valid!");
53 TObjString
tag(
"CovMxManager");
54 dir->WriteTObject(&
tag);
59 const ISyst* syst,
size_t patch) {
64 err <<
"Vector of samples passed to AddSystematic is inconsistent " 65 <<
"with the one passed to the CovMxManager constructor.";
67 throw runtime_error(err.str());
69 for (
size_t i = 0;
i < samples.size(); ++
i) {
71 throw runtime_error(err.str());
88 for (
size_t i = 1;
i <= samples.size(); ++
i) {
92 size_t nPatches = (samples.size() * (samples.size()+1)) / 2;
93 if (patch >= nPatches) {
95 err <<
"Requested patch number " << patch+1 <<
", but there are only" 96 << nPatches <<
" patches in this matrix.";
97 throw runtime_error(err.str());
103 for (
size_t iS = 0; iS < samples.size(); ++iS) {
104 for (
size_t jS = iS; jS < samples.size(); ++jS) {
107 if (counter < patch) {
116 size_t iNBins = samples[iS].GetBinning().NBins();
117 size_t jNBins = samples[jS].GetBinning().NBins();
120 size_t iNComps = iComps.size();
121 size_t jNComps = jComps.size();
124 double iPOT = samples[iS].GetPOT();
125 double jPOT = samples[jS].GetPOT();
128 msg <<
"Producing " << syst->
ShortName() <<
" matrix for samples " 129 << samples[iS].GetName() <<
" and " << samples[jS].GetName();
134 vector< vector<double> > iNom(iNComps, vector<double>(iNBins));
135 vector< vector<double> > jNom(jNComps, vector<double>(jNBins));
138 for (
size_t iC = 0; iC < iNComps; ++iC) {
140 get<1>(iComps[iC]), get<2>(iComps[iC])).
ToTH1(iPOT);
141 for (
size_t iB = 0; iB < iNBins; ++iB) {
142 iNom[iC][iB] = iH->GetBinContent(iB+1);
148 for (
size_t jC = 0; jC < jNComps; ++jC) {
150 get<1>(jComps[jC]), get<2>(jComps[jC])).
ToTH1(jPOT);
151 for (
size_t jB = 0; jB < jNBins; ++jB) {
152 jNom[jC][jB] = jH->GetBinContent(jB+1);
163 for (
size_t u = 0;
u < nUniv; ++
u) {
165 p.SetProgress((
double(
u)+1.)/(
double(nUniv)+1.));
169 double val = rnd.Gaus();
172 cout <<
"Key syst " << tmp->
ShortName() <<
" is one-sided! Throwing positive shift." <<
endl;
178 vector< vector<double> > iShift(iNComps, vector<double>(iNBins));
179 vector< vector<double> > jShift(jNComps, vector<double>(jNBins));
182 for (
size_t iC = 0; iC < iNComps; ++iC) {
184 get<0>(iComps[iC]), get<1>(iComps[iC]), get<2>(iComps[iC])).
ToTH1(iPOT);
185 for (
size_t iB = 0; iB < iNBins; ++iB) {
186 iShift[iC][iB] = iH->GetBinContent(iB+1);
192 for (
size_t jC = 0; jC < jNComps; ++jC) {
194 get<0>(jComps[jC]), get<1>(jComps[jC]), get<2>(jComps[jC])).
ToTH1(jPOT);
195 for (
size_t jB = 0; jB < jNBins; ++jB) {
196 jShift[jC][jB] = jH->GetBinContent(jB+1);
202 vector< vector<double> > iFracShift(iNComps, vector<double>(iNBins));
203 vector< vector<double> > jFracShift(jNComps, vector<double>(jNBins));
206 for (
size_t iC = 0; iC < iNComps; ++iC) {
207 for (
size_t iB = 0; iB < iNBins; ++iB) {
208 if (iNom[iC][iB] > 0) {
209 iFracShift[iC][iB] = ((iShift[iC][iB]+0.01)/(iNom[iC][iB]+0.01)) - 1;
212 iFracShift[iC][iB] = 0;
218 for (
size_t jC = 0; jC < jNComps; ++jC) {
219 for (
size_t jB = 0; jB < jNBins; ++jB) {
220 if (jNom[jC][jB] > 0) {
221 jFracShift[jC][jB] = ((jShift[jC][jB]+0.01)/(jNom[jC][jB]+0.01)) - 1;
224 jFracShift[jC][jB] = 0;
230 for (
size_t iC = 0; iC < iNComps; ++iC) {
231 for (
size_t jC = 0; jC < jNComps; ++jC) {
232 for (
size_t iB = 0; iB < iNBins; ++iB) {
233 for(
size_t jB = 0; jB < jNBins; ++jB) {
234 matrix[iC][jC](iB, jB) += (iFracShift[iC][iB] * jFracShift[jC][jB]) * (1.0 / (
float)nUniv);
266 vector<Sample> samples, vector<const ISyst*>
systs)
269 if (samples.empty())
throw runtime_error(
"Sample vector is empty!");
270 if (systs.empty())
throw runtime_error(
"Systematics vector is empty!");
273 CovarianceMatrix mxTmp(samples);
274 size_t nBinsFull = mxTmp.GetFullBinning().NBins();
275 MatrixXd mxFull(nBinsFull, nBinsFull);
279 cout <<
"Making covariance matrix for systematics:" <<
endl;
280 for (
const ISyst* syst : systs) {
281 cout <<
" " << syst->ShortName() <<
endl;
285 Progress p(
"Producing covariance matrix");
289 int max = systs.size() * samples.size() * samples.size();
290 for (
size_t iSyst = 0; iSyst < systs.size(); ++iSyst) {
292 LoadSyst(systs[iSyst]->ShortName());
295 for (
size_t iS = 0; iS < samples.size(); ++iS) {
296 size_t iNBins = samples[iS].GetBinning().NBins();
298 size_t iNComps = iComps.size();
301 for (
size_t i = 0;
i < iS; ++
i) {
306 for (
size_t jS = 0; jS < samples.size(); ++jS) {
309 size_t jNBins = samples[jS].GetBinning().NBins();
311 size_t jNComps = jComps.size();
314 for (
size_t i = 0;
i < jS; ++
i) {
320 if (mx.size() == 0) {
322 err <<
"Manager not set up for samples " << samples[iS].GetName()
323 <<
" and " << samples[jS].GetName() <<
endl;
324 throw runtime_error(err.str());
327 for (
size_t iC = 0; iC < iNComps; ++iC) {
328 for (
size_t jC = 0; jC < jNComps; ++jC) {
329 for (
size_t iB = 0; iB < iNBins; ++iB) {
330 size_t iBin = iOffset + (iC * iNBins) + iB;
331 for (
size_t jB = 0; jB < jNBins; ++jB) {
332 size_t jBin = jOffset + (jC * jNBins) + jB;
333 mxFull(iBin,jBin) += mx[iC][jC](iB,jB);
344 return make_unique<CovarianceMatrix>(mxFull, samples);
364 TDirectory* systsDir =
fTopDir->GetDirectory(
"systs");
365 if (!systsDir) systsDir =
fTopDir->mkdir(
"systs");
366 TDirectory* dir = systsDir->GetDirectory(mxName.c_str());
367 if (!dir) dir = systsDir->mkdir(mxName.c_str());
372 ostringstream samples;
373 samples << mx.first.first.GetID() <<
"_" << mx.first.second.GetID();
376 TDirectory* patchDir = dir->GetDirectory(samples.str().c_str());
377 if (!patchDir) patchDir = dir->mkdir(samples.str().c_str());
381 int b1 = mx.first.first.GetBinning().NBins();
383 int b2 = mx.first.second.GetBinning().NBins();
384 vector<int>
bins = {
c1, b1,
c2, b2 };
385 patchDir->WriteObject(&bins,
"bins");
387 for (
int iC = 0; iC <
c1; ++iC) {
388 for (
int jC = 0; jC <
c2; ++jC) {
391 key << iC <<
"_" << jC;
392 patchDir->WriteTObject(&mx.second[iC][jC], key.str().c_str());
410 throw runtime_error(
"CovMxManager directory not valid!");
412 TDirectory* dir =
fTopDir->GetDirectory(
"systs")->GetDirectory(mxName.c_str());
415 err <<
"CovMxManager directory for systematic " << mxName
417 throw runtime_error(err.str());
419 TIter
next(dir->GetListOfKeys());
421 while ((key = (TKey*)
next())) {
423 TDirectory* patchDir = (TDirectory*)key->ReadObj();
424 string name(key->GetName());
428 size_t end =
name.find(
"_", start);
429 int id = stoi(
name.substr(start, end));
433 id = stoi(
name.substr(start));
440 patchDir->GetObject(
"bins", bins);
441 int c1 = bins->at(0);
442 int b1 = bins->at(1);
443 int c2 = bins->at(2);
444 int b2 = bins->at(3);
450 for (
int iC = 0; iC <
c1; ++iC) {
451 for (
int jC = 0; jC <
c2; ++jC) {
452 ostringstream mxname;
453 mxname << iC <<
"_" << jC;
455 if (!mxtemp)
throw runtime_error(
"Failed to load matrix with name "+mxname.str());
456 matrix[iC][jC] = *mxtemp;
470 Patch matrixTrans(c2, vector< TMatrixD >(c1,
TMatrixD(b2, b1)));
471 for (
int iC = 0; iC <
c1; ++iC) {
472 for (
int jC = 0; jC <
c2; ++jC) {
473 matrixTrans[jC][iC].Transpose(matrix[iC][jC]);
476 fCovMx[pairTrans] = matrixTrans;
std::unique_ptr< covmx::CovarianceMatrix > GetCovarianceMatrix(std::vector< covmx::Sample > samples, std::vector< const ISyst * > systs)
Cuts and Vars for the 2020 FD DiF Study.
vector< Component > GetComponents(Sample sample)
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
fvar< T > fabs(const fvar< T > &x)
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.
virtual const std::string & ShortName() const final
The name printed out to the screen.
Simple record of shifts applied to systematic parameters.
void LoadSyst(std::string mxName)
Function to load matrix for a given systematic.
virtual Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Version of OscCalcSterile that always returns probability of 1.
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Encapsulate code to systematically shift a caf::SRProxy.
std::pair< Sample, Sample > SamplePair
~CovMxManager()
Destructor.
std::vector< unsigned int > fSampleIDs
TMatrixT< double > TMatrixD
std::map< covmx::SamplePair, covmx::Patch > fCovMx
void AddSystematic(std::vector< covmx::Sample > &samples, const ISyst *syst, size_t patch)
void SetProgress(double frac)
Update the progress fraction between zero and one.
static std::unique_ptr< CovMxManager > LoadFrom(TDirectory *dir)
Standard implementation.
A simple ascii-art progress bar.
Standard interface to all prediction techniques.
std::vector< std::vector< TMatrixD > > Patch
Prevent histograms being added to the current directory.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
void Done()
Call this when action is completed.
void SetShift(const ISyst *syst, double shift, bool force=false)
CovMxManager(TDirectory *dir, std::vector< covmx::Sample > samples)
void SaveSyst(std::string mxName)
Function to save matrix for a given systematic.