MatrixDeterminant.C
Go to the documentation of this file.
1 /// MatrixDeterminant.C by J. Hewes (jhewes15@fnal.gov)
2 /// Test how the fractional matrix determinant changes in parameter space.
3 
4 #include "CAFAna/Core/Progress.h"
9 
10 #include <cmath>
11 
12 using namespace ana;
13 
14 TVectorD hist2vec(std::vector<TH1D*> h) {
15  size_t nElements = 0;
16  for (TH1D* hTemp : h) nElements += hTemp->GetNbinsX();
17  TVectorD v(nElements);
18  size_t count = 0;
19  for (TH1D* hTemp : h) {
20  for (size_t i = 0; i < hTemp->GetNbinsX(); ++i) {
21  v(count) = hTemp->GetBinContent(i+1);
22  ++count;
23  }
24  }
25  for (TH1D* hTemp : h) delete hTemp;
26  return v;
27 }
28 
30 
31  DontAddDirectory guard;
32 
33  // Define exposures
34  double potND = kAna2018SensitivityFHCNDPOT;
35  double potFD = kAna2018FHCPOT;
36 
37  auto calc = DefaultSterileCalc(4);
39  IPrediction* predND = LoadPrediction(kNusFHCNearDet, "nosysts", kUseCVMFS);
40  IPrediction* predFD = LoadPrediction(kNusFHCFarDet, "nosysts", kUseCVMFS);
41  CovarianceMatrix* mx = GetFullMatrix();
42 
43  Binning xbins = Binning::LogUniform(50, 1e-4, 1);
44  Binning ybins = Binning::LogUniform(50, 1e-2, 100);
45 
46  TH2D* hist = new TH2D("h", "h", xbins.NBins(), &xbins.Edges()[0],
47  ybins.NBins(), &ybins.Edges()[0]);
48  Progress p("Calculating determinants");
49  for (size_t idx = 1; idx <= hist->GetNbinsX(); ++idx) {
50  for (size_t idy = 1; idy <= hist->GetNbinsY(); ++idy) {
51  kFitSinSqTheta24Sterile.SetValue(calc, hist->GetXaxis()->GetBinCenter(idx));
52  kFitDmSq41Sterile.SetValue(calc, hist->GetYaxis()->GetBinCenter(idy));
53  std::vector<TH1D*> hists = {predND->Predict(calc).ToTH1(potND), predFD->Predict(calc).ToTH1(potFD)};
54  TVectorD vec = hist2vec(hists);
55  mx->PredictCovMx({predND,predFD}, {potND,potFD}, calc);
56  TMatrixD tmx = mx->GetCovMxRelative(vec);
57  double det = tmx.Determinant();
58  // std::cout << "Determinant in bin " << idx << ", " << idy << " is " << det << std::endl;
59  hist->SetBinContent(idx, idy, log(det));
60  p.SetProgress(((50*(idx-1))+idy)/2500.);
61  } // for bin y
62  } // for bin x
63  p.Done();
64 
65  TFile* fOut = TFile::Open("matrix_determinant.root", "recreate");
66  fOut->WriteTObject(hist, "determinant");
67  fOut->WriteTObject(mx->GetFullCovMxTH2(), "fullmx");
68 
69 }
70 
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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
const char * p
Definition: xmltok.h:285
void SetNus20Params(osc::IOscCalcSterile *calc, std::string type)
Definition: Calcs.cxx:29
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
TString hists[nhists]
Definition: bdt_com.C:3
osc::OscCalcDumb calc
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
const int xbins
Definition: MakePlots.C:82
void MatrixDeterminant()
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
const std::vector< double > & Edges() const
Definition: Binning.h:34
const int ybins
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
static Binning LogUniform(int n, double lo, double hi)
Definition: Binning.cxx:118
int NBins() const
Definition: Binning.h:29
TVectorD hist2vec(std::vector< TH1D * > h)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
A simple ascii-art progress bar.
Definition: Progress.h:9
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
const FitDmSq41Sterile kFitDmSq41Sterile