CovMxManager.cxx
Go to the documentation of this file.
1 // CAFAna includes
3 #include "CAFAna/Core/Progress.h"
8 
10 
11 // ROOT includes
12 #include "TKey.h"
13 #include "TMatrixD.h"
14 
15 #include <Eigen/Dense>
16 
17 // C++ includes
18 #include <iostream>
19 #include <numeric>
20 
21 namespace ana {
22 
23  using std::back_inserter;
24  using std::begin;
25  using std::cout;
26  using std::end;
27  using std::endl;
28  using std::get;
29  using std::make_pair;
30  using std::make_unique;
31  using std::ostringstream;
32  using std::runtime_error;
33  using std::stoi;
34  using std::string;
35  using std::transform;
36  using std::unique_ptr;
37  using std::vector;
38 
39  using Eigen::MatrixXd;
40 
41  using covmx::CovarianceMatrix;
42  using covmx::Patch;
43  using covmx::Sample;
44  using covmx::SamplePair;
45 
46  //----------------------------------------------------------------------
47  CovMxManager::CovMxManager(TDirectory* dir, std::vector<Sample> samples)
48  : fTopDir(dir) {
49 
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);
55  }
56 
57  //----------------------------------------------------------------------
58  void CovMxManager::AddSystematic(vector<Sample>& samples,
59  const ISyst* syst, size_t patch) {
60 
61  // Check samples are consistent
62  if (true) { // Let objects in this code block fall out of scope
63  ostringstream err;
64  err << "Vector of samples passed to AddSystematic is inconsistent "
65  << "with the one passed to the CovMxManager constructor.";
66  if (samples.size() != fSampleIDs.size()) {
67  throw runtime_error(err.str());
68  }
69  for (size_t i = 0; i < samples.size(); ++i) {
70  if (samples[i].GetID() != fSampleIDs[i]) {
71  throw runtime_error(err.str());
72  }
73  } // for sample i
74  } // if true
75 
76  // Trivial oscillation calculator
78 
79  // Making a matrix across many samples can be computationally intensive,
80  // so we want to be able to break it up into smaller chunks. However,
81  // because it's a symmetric matrix and we only actually generate in a
82  // triangular manner and then mirror when necessary, splitting it up
83  // evenly isn't entirely trivial.
84 
85  // Overall we will be generating n factorial matrices, where n is the
86  // number of samples. That's what the "total" variable is
87  size_t total = 0;
88  for (size_t i = 1; i <= samples.size(); ++i) {
89  total += i;
90  }
91 
92  size_t nPatches = (samples.size() * (samples.size()+1)) / 2;
93  if (patch >= nPatches) {
94  ostringstream err;
95  err << "Requested patch number " << patch+1 << ", but there are only"
96  << nPatches << " patches in this matrix.";
97  throw runtime_error(err.str());
98  }
99 
100  // Loop over samples
101  size_t counter = 0;
102 
103  for (size_t iS = 0; iS < samples.size(); ++iS) {
104  for (size_t jS = iS; jS < samples.size(); ++jS) {
105 
106  // Skip matrices until we get to the first one we want
107  if (counter < patch) {
108  ++counter;
109  continue;
110  }
111 
112  IPrediction* iP = samples[iS].GetPrediction();
113  IPrediction* jP = samples[jS].GetPrediction();
114 
115  // Number of bins & components for each of the current samples
116  size_t iNBins = samples[iS].GetBinning().NBins();
117  size_t jNBins = samples[jS].GetBinning().NBins();
118  auto iComps = covmx::GetComponents(samples[iS]);
119  auto jComps = covmx::GetComponents(samples[jS]);
120  size_t iNComps = iComps.size();
121  size_t jNComps = jComps.size();
122 
123  // Get exposures
124  double iPOT = samples[iS].GetPOT();
125  double jPOT = samples[jS].GetPOT();
126 
127  ostringstream msg;
128  msg << "Producing " << syst->ShortName() << " matrix for samples "
129  << samples[iS].GetName() << " and " << samples[jS].GetName();
130  Progress p(msg.str());
131  size_t nUniv = 1e5;
132 
133  // Nominal spectra
134  vector< vector<double> > iNom(iNComps, vector<double>(iNBins));
135  vector< vector<double> > jNom(jNComps, vector<double>(jNBins));
136 
137  // Populate nominal spectra for sample i
138  for (size_t iC = 0; iC < iNComps; ++iC) {
139  TH1D* iH = iP->PredictComponent(calc, get<0>(iComps[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);
143  }
144  delete iH;
145  } // for component i
146 
147  // Populate nominal spectra for sample j
148  for (size_t jC = 0; jC < jNComps; ++jC) {
149  TH1D* jH = jP->PredictComponent(calc, get<0>(jComps[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);
153  }
154  delete jH;
155  } // for component i
156 
157  // Prepare key and matrix vector inputs
158  SamplePair key = make_pair(samples[iS], samples[jS]);
159  Patch matrix(iNComps, vector< TMatrixD >(jNComps, TMatrixD(iNBins, jNBins)));
160 
161  TRandom3 rnd(0);
162 
163  for (size_t u = 0; u < nUniv; ++u) {
164 
165  p.SetProgress((double(u)+1.)/(double(nUniv)+1.));
166 
167  // Throw the systematic
168  SystShifts shift;
169  double val = rnd.Gaus();
170  const KeySyst* tmp = dynamic_cast<const KeySyst*>(syst);
171  if (tmp && tmp->IsOneSided()) {
172  cout << "Key syst " << tmp->ShortName() << " is one-sided! Throwing positive shift." << endl;
173  val = fabs(val);
174  }
175  shift.SetShift(syst, val);
176 
177  // Get systematic shifts
178  vector< vector<double> > iShift(iNComps, vector<double>(iNBins));
179  vector< vector<double> > jShift(jNComps, vector<double>(jNBins));
180 
181  // Populate shifted spectra for sample i
182  for (size_t iC = 0; iC < iNComps; ++iC) {
183  TH1D* iH = iP->PredictComponentSyst(calc, samples[iS].GetSystShifts(shift),
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);
187  }
188  delete iH;
189  } // for component i
190 
191  // Populate shifted spectra for sample j
192  for (size_t jC = 0; jC < jNComps; ++jC) {
193  TH1D* jH = jP->PredictComponentSyst(calc, samples[jS].GetSystShifts(shift),
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);
197  }
198  delete jH;
199  } // for component j
200 
201  // Get fractional shifts
202  vector< vector<double> > iFracShift(iNComps, vector<double>(iNBins));
203  vector< vector<double> > jFracShift(jNComps, vector<double>(jNBins));
204 
205  // Populate fractional shift for sample i
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;
210  }
211  else {
212  iFracShift[iC][iB] = 0;
213  }
214  } // for bin i
215  } // for component i
216 
217  // Populate fractional shift for sample j
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;
222  }
223  else {
224  jFracShift[jC][jB] = 0;
225  }
226  } // for bin j
227  } // for component j
228 
229  // Now populate matrix
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);
235  } // for bin j
236  } // for bin i
237  } // for component j
238  } // for component i
239  } // for universe u
240 
241  p.Done();
242  fCovMx.insert(make_pair(key, matrix));
243 
244  // Save the matrices and remove them from memory
245  SaveSyst(syst->ShortName());
246  fCovMx.clear();
247 
248  return; // Done processing this patch
249 
250  } // for sample j
251  } // for sample i
252 
253  } // function CovMxManager::AddSystematic
254 
255  //----------------------------------------------------------------------
256  /// Constructor for LoadFrom implementation
257  CovMxManager::CovMxManager(TDirectory* dir)
258  : fTopDir(dir) {}
259 
260  //----------------------------------------------------------------------
261  /// Destructor
263 
264  //----------------------------------------------------------------------
265  unique_ptr<CovarianceMatrix> CovMxManager::GetCovarianceMatrix(
266  vector<Sample> samples, vector<const ISyst*> systs)
267  {
268  // Quick sanity check
269  if (samples.empty()) throw runtime_error("Sample vector is empty!");
270  if (systs.empty()) throw runtime_error("Systematics vector is empty!");
271 
272  // Temporary matrix just to figure out binning
273  CovarianceMatrix mxTmp(samples);
274  size_t nBinsFull = mxTmp.GetFullBinning().NBins();
275  MatrixXd mxFull(nBinsFull, nBinsFull);
276  mxFull.setZero();
277 
278  cout << endl;
279  cout << "Making covariance matrix for systematics:" << endl;
280  for (const ISyst* syst : systs) {
281  cout << " " << syst->ShortName() << endl;
282  }
283  cout << endl;
284 
285  Progress p("Producing covariance matrix");
286 
287  // Loop over each systematic
288  int counter = 0;
289  int max = systs.size() * samples.size() * samples.size();
290  for (size_t iSyst = 0; iSyst < systs.size(); ++iSyst) {
291 
292  LoadSyst(systs[iSyst]->ShortName());
293 
294  // Loop over samples in i
295  for (size_t iS = 0; iS < samples.size(); ++iS) {
296  size_t iNBins = samples[iS].GetBinning().NBins();
297  auto iComps = covmx::GetComponents(samples[iS]);
298  size_t iNComps = iComps.size();
299  // Figure out offset in i
300  size_t iOffset(0);
301  for (size_t i = 0; i < iS; ++i) {
302  // Offset the full matrix index by the total number of matrix elements in preceding samples
303  iOffset += samples[i].GetBinning().NBins() * covmx::GetComponents(samples[i]).size();
304  }
305  // Loop over samples in j
306  for (size_t jS = 0; jS < samples.size(); ++jS) {
307  p.SetProgress((double)counter/(double)max);
308  ++counter;
309  size_t jNBins = samples[jS].GetBinning().NBins();
310  auto jComps = covmx::GetComponents(samples[jS]);
311  size_t jNComps = jComps.size();
312  // Figure out offset in j
313  size_t jOffset(0);
314  for (size_t i = 0; i < jS; ++i) {
315  // Offset the full matrix index by the total number of matrix elements in preceding samples
316  jOffset += samples[i].GetBinning().NBins() * covmx::GetComponents(samples[i]).size();
317  }
318  SamplePair key = make_pair(samples[iS], samples[jS]);
319  Patch mx = fCovMx[key];
320  if (mx.size() == 0) {
321  ostringstream err;
322  err << "Manager not set up for samples " << samples[iS].GetName()
323  << " and " << samples[jS].GetName() << endl;
324  throw runtime_error(err.str());
325  }
326  // Loop over individual components and bins within samples
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);
334  } // for bin j
335  } // for bin i
336  } // for component j
337  } // for component i
338  } // for sample j
339  } // for sample i
340  } // for syst
341 
342  p.Done();
343 
344  return make_unique<CovarianceMatrix>(mxFull, samples);
345 
346  } // function CovMxManager::GetCovarianceMatrix
347 
348  //----------------------------------------------------------------------
349  /// Standard implementation
350  unique_ptr<CovMxManager> CovMxManager::LoadFrom(
351  TDirectory* dir)
352  {
353  unique_ptr<CovMxManager> ret(new CovMxManager(dir));
354 
355  return move(ret);
356 
357  } // function CovMxManager::LoadFrom
358 
359  //----------------------------------------------------------------------
360  /// Function to save matrix for a given systematic
361  void CovMxManager::SaveSyst(string mxName)
362  {
363  // Create directory for this systematic
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());
368 
369  for (auto mx : fCovMx) {
370 
371  // Write a key for this pair of samples
372  ostringstream samples;
373  samples << mx.first.first.GetID() << "_" << mx.first.second.GetID();
374 
375  // Create a directory
376  TDirectory* patchDir = dir->GetDirectory(samples.str().c_str());
377  if (!patchDir) patchDir = dir->mkdir(samples.str().c_str());
378 
379  // Save number of components & bins for each sample pair
380  int c1 = covmx::GetComponents(mx.first.first).size();
381  int b1 = mx.first.first.GetBinning().NBins();
382  int c2 = covmx::GetComponents(mx.first.second).size();
383  int b2 = mx.first.second.GetBinning().NBins();
384  vector<int> bins = { c1, b1, c2, b2 };
385  patchDir->WriteObject(&bins, "bins");
386  // Save matrix patch for each component
387  for (int iC = 0; iC < c1; ++iC) {
388  for (int jC = 0; jC < c2; ++jC) {
389  // Save each matrix with the corresponding key
390  ostringstream key;
391  key << iC << "_" << jC;
392  patchDir->WriteTObject(&mx.second[iC][jC], key.str().c_str());
393  } // for component j
394  } // for component i
395  } // for sample pair mx
396 
397  } // function CovMxManager::SaveSyst
398 
399  //----------------------------------------------------------------------
400  /// Function to load matrix for a given systematic
402  {
403  DontAddDirectory guard;
404  fCovMx.clear();
405 
406  // We're loading the covariance matrix for a given systematic.
407  // This means looping over all the samples in this systematic's
408  // directory, and loading the corresponding matrix patch.
409  if (!fTopDir) {
410  throw runtime_error("CovMxManager directory not valid!");
411  }
412  TDirectory* dir = fTopDir->GetDirectory("systs")->GetDirectory(mxName.c_str());
413  if (!dir) {
414  ostringstream err;
415  err << "CovMxManager directory for systematic " << mxName
416  << " not valid!";
417  throw runtime_error(err.str());
418  }
419  TIter next(dir->GetListOfKeys());
420  TKey* key;
421  while ((key = (TKey*)next())) {
422 
423  TDirectory* patchDir = (TDirectory*)key->ReadObj();
424  string name(key->GetName());
425 
426  // Parse key string to get sample information
427  size_t start = 0;
428  size_t end = name.find("_", start);
429  int id = stoi(name.substr(start, end));
430  Sample s1(id);
431 
432  start = end + 1;
433  id = stoi(name.substr(start));
434  Sample s2(id);
435 
436  // Now we make a sample pair, and create the patch of covariance matrix
437  // corresponding to this pair of samples
438  SamplePair pair = make_pair(s1, s2);
439  vector<int>* bins;
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);
445 
446  Patch matrix(c1, vector< TMatrixD >(c2, TMatrixD(b1, b2)));
447 
448  // Now we loop over and load all the matrices for the various beam
449  // components which make up this patch of covariance matrix
450  for (int iC = 0; iC < c1; ++iC) {
451  for (int jC = 0; jC < c2; ++jC) {
452  ostringstream mxname;
453  mxname << iC << "_" << jC;
454  TMatrixD* mxtemp = (TMatrixD*)patchDir->Get(mxname.str().c_str());
455  if (!mxtemp) throw runtime_error("Failed to load matrix with name "+mxname.str());
456  matrix[iC][jC] = *mxtemp;
457  delete mxtemp;
458  } // for component j
459  } // for component i
460 
461  // Add to map
462  fCovMx[pair] = matrix;
463 
464  // If we're on an off-diagonal patch, then there's a corresponding
465  // mirrored patch on the other side of the diagonal. It's exactly
466  // equivalent to this one, so rather than storing redundant information
467  // we just mirror this patch into that one.
468  if (s1 != s2) {
469  SamplePair pairTrans = make_pair(s2, s1);
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]);
474  } // for component j
475  } // for component i
476  fCovMx[pairTrans] = matrixTrans;
477  } // if transposing
478 
479  delete patchDir;
480 
481  } // for sample pair key in directory
482 
483  delete dir;
484 
485  } // function CovMxManager::LoadSyst
486 
487 } // namespace ana
488 
const XML_Char * name
Definition: expat.h:151
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.
Definition: vars.h:6
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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:209
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void LoadSyst(std::string mxName)
Function to load matrix for a given systematic.
const char * p
Definition: xmltok.h:285
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.
Float_t tmp
Definition: plot.C:36
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
TDirectory * fTopDir
Definition: CovMxManager.h:55
std::pair< Sample, Sample > SamplePair
Definition: CovMxManager.h:24
c2
Definition: demo5.py:33
const XML_Char * s
Definition: expat.h:262
bool IsOneSided() const
Definition: CovMxSysts.h:41
~CovMxManager()
Destructor.
std::vector< unsigned int > fSampleIDs
Definition: CovMxManager.h:58
TMatrixT< double > TMatrixD
Definition: Utilities.h:16
std::map< covmx::SamplePair, covmx::Patch > fCovMx
Definition: CovMxManager.h:57
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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.
Definition: Progress.cxx:41
TDirectory * dir
Definition: macro.C:5
static std::unique_ptr< CovMxManager > LoadFrom(TDirectory *dir)
Standard implementation.
A simple ascii-art progress bar.
Definition: Progress.h:9
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::vector< std::vector< TMatrixD > > Patch
Definition: CovMxManager.h:23
Prevent histograms being added to the current directory.
Definition: Utilities.h:62
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void Done()
Call this when action is completed.
Definition: Progress.cxx:91
void next()
Definition: show_event.C:84
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
CovMxManager(TDirectory *dir, std::vector< covmx::Sample > samples)
void SaveSyst(std::string mxName)
Function to save matrix for a given systematic.