ScaleCovarianceMatrix.C
Go to the documentation of this file.
1 // script to pull the prod4 covariance matrix and svae it in a file
2 // note that this needs quite a few local modifications (or to be run
3 // in the old prod4 branch)
4 //
5 // author: Adam Lister
6 // date : 2020-04-09
7 // email : adam.lister@wisc.edu
8 
10 #include "CAFAna/Core/Utilities.h"
11 #include "CAFAna/Core/Sample.h"
12 #include "CAFAna/Analysis/Calcs.h"
18 
24 
25 #include "OscLib/OscCalcSterile.h"
26 
27 std::vector<double> hist2vec(std::vector<TH1D*> h) {
28  size_t nElements = 0;
29  for (TH1D* hTemp : h) nElements += hTemp->GetNbinsX();
30  std::vector<double> v(nElements);
31  size_t count = 0;
32  for (TH1D* hTemp : h) {
33  for (size_t i = 0; i < hTemp->GetNbinsX(); ++i) {
34  std::cout << i << "/" << hTemp->GetNbinsX() << std::endl;
35  v.at(count) = hTemp->GetBinContent(i+1);
36  ++count;
37  }
38  }
39  for (TH1D* hTemp : h) delete hTemp;
40  return v;
41 }
42 
44 
45  // setup oscillation calculator
46  auto calc = DefaultSterileCalc(4);
47 
48  SetNus20Params(calc); // initalise to null hypothesis (no sterile oscillations)
50 
51  // get samples (i.e nus fhc nd, nus fhc fardet)
52  std::vector<covmx::Sample> samples = GetSamplesFromOptString("nus_fhc_neardet_fardet");
53 
54  // loop samples, get each prediction object, and push back TH1
55  // prediction objects to a spectra vector
56  std::vector<TH1D*> spectras;
57  for (covmx::Sample& sample : samples) {
58  IPrediction* pred = LoadPrediction(sample, "nosysts", false);
59  sample.SetPrediction(pred);
60  spectras.push_back(pred->Predict(calc).ToTH1(sample.GetPOT()));
61  }
62 
63  // get a vector from the full spectra
64  std::vector<double> vec = hist2vec(spectras);
65 
66  // this is a vector of CovarianceMatrix pointers
67  auto matrices = GetMatrices(samples, false);
68 
69  // get the full covariance matrix and
70  // fractional covariance matrix
71  std::vector<TMatrixD> totcovmx;
72  std::vector<TMatrixD> fraccovmx;
73  for (covmx::CovarianceMatrix* mx : matrices) {
74  totcovmx .push_back(mx->Predict(samples, calc));
75  fraccovmx.push_back(mx->GetCovMxRelative(vec));
76  }
77 
78  TFile* covmxfile = new TFile("covmx.root", "recreate");
79  covmxfile->cd();
80  TMatrixD fulltotcovmx (totcovmx.at(0).GetNrows(),
81  totcovmx.at(0).GetNcols());
82  TMatrixD fullfraccovmx(totcovmx.at(0).GetNrows(),
83  totcovmx.at(0).GetNcols());
84 
85  for (size_t i = 0; i < totcovmx.size(); ++i){
86  fulltotcovmx = fulltotcovmx + totcovmx.at(i);
87  fullfraccovmx = fullfraccovmx + fraccovmx.at(i);
88  }
89  fulltotcovmx .Write("covmx");
90  fullfraccovmx.Write("fraccovmx");
91 }
void ScaleCovarianceMatrix()
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
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
std::vector< double > hist2vec(std::vector< TH1D * > h)
void PrintOscParams(osc::OscCalcSterile *calc)
osc::OscCalcDumb calc
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
Eigen::VectorXd vec
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
OStream cout
Definition: OStream.cxx:6
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
std::vector< covmx::Sample > GetSamplesFromOptString(TString optString)
Function to take an option TString and return a vector of associated covmx::Samples.
Definition: Utilities.h:384
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...