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"
17 
23 
24 std::vector<double> hist2vec(std::vector<TH1D*> h) {
25  size_t nElements = 0;
26  for (TH1D* hTemp : h) nElements += hTemp->GetNbinsX();
27  std::vector<double> v(nElements);
28  size_t count = 0;
29  for (TH1D* hTemp : h) {
30  for (size_t i = 0; i < hTemp->GetNbinsX(); ++i) {
31  std::cout << i << "/" << hTemp->GetNbinsX() << std::endl;
32  v.at(count) = hTemp->GetBinContent(i+1);
33  ++count;
34  }
35  }
36  for (TH1D* hTemp : h) delete hTemp;
37  return v;
38 }
39 
41 
42  // setup oscillation calculator
43  auto calc = DefaultSterileCalc(4);
44 
45  SetNus20Params(calc); // initalise to null hypothesis (no sterile oscillations)
47 
48  // get samples (i.e nus fhc nd, nus fhc fardet)
49  std::vector<covmx::Sample> samples = GetSamplesFromOptString("nus_fhc_neardet_fardet");
50 
51  // loop samples, get each prediction object, and push back TH1
52  // prediction objects to a spectra vector
53  std::vector<TH1D*> spectras;
54  for (covmx::Sample& sample : samples) {
55  IPrediction* pred = LoadPrediction(sample, "nosysts", false);
56  sample.SetPrediction(pred);
57  spectras.push_back(pred->Predict(calc).ToTH1(sample.GetPOT()));
58  }
59 
60  // get a vector from the full spectra
61  std::vector<double> vec = hist2vec(spectras);
62 
63  // this is a vector of CovarianceMatrix pointers
64  auto matrices = GetMatrices(samples, false);
65 
66  // get the full covariance matrix and
67  // fractional covariance matrix
68  std::vector<TMatrixD> totcovmx;
69  std::vector<TMatrixD> fraccovmx;
70  for (covmx::CovarianceMatrix* mx : matrices) {
71  totcovmx .push_back(mx->Predict(samples, calc));
72  fraccovmx.push_back(mx->GetCovMxRelative(vec));
73  }
74 
75  TFile* covmxfile = new TFile("covmx.root", "recreate");
76  covmxfile->cd();
77  TMatrixD fulltotcovmx (totcovmx.at(0).GetNrows(),
78  totcovmx.at(0).GetNcols());
79  TMatrixD fullfraccovmx(totcovmx.at(0).GetNrows(),
80  totcovmx.at(0).GetNcols());
81 
82  for (size_t i = 0; i < totcovmx.size(); ++i){
83  fulltotcovmx = fulltotcovmx + totcovmx.at(i);
84  fullfraccovmx = fullfraccovmx + fraccovmx.at(i);
85  }
86  fulltotcovmx .Write("covmx");
87  fullfraccovmx.Write("fraccovmx");
88 }
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:176
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
std::vector< double > hist2vec(std::vector< TH1D * > h)
void SetNus20Params(osc::IOscCalcSterile *calc, std::string type)
Definition: Calcs.cxx:29
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
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:279
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 ...