DecorrelateFD.C
Go to the documentation of this file.
1 /// Macro to decorrelate far detector errors against near detector
2 
5 
6 using namespace ana;
7 
9 {
10  DontAddDirectory guard;
11 
12  // Define exposures
13  double potND = kAna2018SensitivityFHCNDPOT;
14  double potFD = kAna2018FHCPOT;
15 
16  // Load covariance matrix and predictions
17  CovarianceMatrix mx({kNCNDBinning,kNCFDBinning});
18  for (auto syst : GetMatrices({kNusFHCNearDet,kNusFHCFarDet}, kUseCVMFS)) mx.AddMatrix(syst);
19  std::vector<IPrediction*> preds = { LoadPrediction(kNusFHCNearDet, "nosysts", kUseCVMFS),
20  LoadPrediction(kNusFHCFarDet, "nosysts", kUseCVMFS) };
21  auto calc = DefaultSterileCalc(4);
23  TMatrixD matrix = mx.PredictCovMx(preds, {potND, potFD}, calc);
24 
25  size_t binningND = kNCNDBinning.NBins();
26  size_t binningFD = kNCFDBinning.NBins();
27 
28  std::vector<double> vFD(binningFD);
29  TH1D* hFD = preds[1]->Predict(calc).ToTH1(potFD);
30  for (size_t i = 0; i < binningND; ++i)
31  vFD[i] = hFD->GetBinContent(i+1);
32 
33  // Split the matrix up into ND, FD and ND-FD quadrants
34  TMatrixD s11(binningFD, binningFD);
35  TMatrixD s12(binningND, binningFD);
36  TMatrixD s21(binningFD, binningND);
37  TMatrixD s22(binningND, binningND);
38 
39  for (size_t i = 0; i < binningND; ++i) {
40  for (size_t j = 0; j < binningND; ++j) {
41  s22(i,j) = matrix(i,j);
42  }
43  }
44 
45  for (size_t i = 0; i < binningND; ++i) {
46  for (size_t j = 0; j < binningFD; ++j) {
47  s12(i,j) = matrix(i,binningND+j);
48  s21(j,i) = matrix(i,binningND+j);
49  }
50  }
51 
52  for (size_t i = 0; i < binningFD; ++i) {
53  for (size_t j = 0; j < binningFD; ++j) {
54  s11(i,j) = matrix(binningND+i,binningND+j);
55  }
56  }
57 
58  std::vector<double> before;
59  for (size_t i = 0; i < binningFD; ++i) {
60  before.push_back(sqrt(matrix(binningND+i,binningND+i))/vFD[i]);
61  }
62 
63  // Carry out the decorrelation procedure
64  TMatrixD tmp(binningFD, binningFD);
65  s22.Invert();
66  s21 *= s22;
67  tmp.Mult(s21, s12);
68  s11 -= tmp;
69 
70  std::vector<double> after;
71  for (size_t i = 0; i < binningFD; ++i) {
72  after.push_back(sqrt(s11(i,i))/vFD[i]);
73  }
74 
75  // Print the reduction in uncertainty
76  for (size_t i = 0; i < binningFD; ++i) {
77  std::cout << "Uncertainty shrinks from " << before[i] << " to " << after[i] << std::endl;
78  }
79 
80 
81  // Draw a histogram of errors before and after
82  TH1D* hBefore = new TH1D("hBefore", ";Energy deposited in FD (GeV);Systematic uncertainty (%)",
84  TH1D* hAfter = new TH1D("hAfter", ";Energy deposited in FD (GeV);Systematic uncertainty (%)",
86  for (size_t i = 0; i < binningFD; ++i) {
87  hBefore->SetBinContent(i+1, 100*before[i]);
88  hAfter->SetBinContent(i+1, 100*after[i]);
89  }
90 
91  // Format and save the plot
92  TCanvas* c = new TCanvas("c", "", 1600, 900);
93  hBefore->SetLineColor(kRed+2);
94  hBefore->SetLineWidth(3);
95  hAfter->SetLineColor(kGreen+2);
96  hAfter->SetLineWidth(3);
97  hBefore->Draw("hist");
98  hAfter->Draw("hist same");
99  c->SaveAs("decorrelated.pdf");
100 
101 }
102 
def after(args)
Definition: fabricate.py:1416
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
T sqrt(T number)
Definition: d0nt_math.hpp:156
const bool kUseCVMFS
Float_t tmp
Definition: plot.C:36
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
osc::OscCalcDumb calc
void DecorrelateFD()
Definition: DecorrelateFD.C:8
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
const double j
Definition: BetheBloch.cxx:29
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
const std::vector< double > & Edges() const
Definition: Binning.h:34
OStream cout
Definition: OStream.cxx:6
const Binning kNCNDBinning
NC custom binning.
Definition: Binning.cxx:104
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
int NBins() const
Definition: Binning.h:29
const Binning kNCFDBinning
Definition: Binning.cxx:105
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)