TestCovMxNew.C
Go to the documentation of this file.
1 /// MakeCovMx.C by J. Hewes (jhewes15@fnal.gov)
2 // Macro to make covariance matrices
3 
4 #include "TCanvas.h"
5 
9 
10 using namespace ana;
11 
12 const std::vector<std::string> ndComps = { "cc_nue_to_nue",
13  "cc_nuebar_to_nuebar",
14  "cc_numu_to_numu",
15  "cc_numubar_to_numubar",
16  "nc" };
17 
18 const std::vector<std::string> fdComps = { "cc_nue_to_nue",
19  "cc_nuebar_to_nuebar",
20  "cc_nue_to_numu",
21  "cc_nuebar_to_numubar",
22  "cc_nue_to_nutau",
23  "cc_nuebar_to_nutaubar",
24  "cc_numu_to_nue",
25  "cc_numubar_to_nuebar",
26  "cc_numu_to_numu",
27  "cc_numubar_to_numubar",
28  "cc_numu_to_nutau",
29  "cc_numubar_to_nutaubar",
30  "nc" };
31 
32 std::vector<TH1D*> GetSpectraFromDir(TDirectory* dir, std::vector<std::string> comps, double pot);
33 
35 
36  try {
37 
38  std::vector<covmx::Sample> samples = GetSamplesFromOptString(opt);
39 
40  for (covmx::Sample& sample : samples) {
41  SetBinning(sample);
42  if (sample.detector == covmx::kNearDet) {
43  sample.SetPOT(kAna2018SensitivityFHCNDPOT);
44  sample.SetLivetime(-1);
45  }
46  else {
47  sample.SetPOT(kAna2018FHCPOT);
48  sample.SetLivetime(kAna2018FHCLivetime);
49  }
50  SetPrediction(sample, false);
51  }
52 
53  // Load systematics
54  std::vector<std::string> systs;
55  for (const ISyst* syst : getAllXsecSysts_2018_RPAFix()) {
56  systs.push_back(syst->ShortName());
57  }
58  for (const ISyst* syst : LoadSysts(samples[0], kUseCVMFS)) {
59  systs.push_back(dynamic_cast<const NuISyst*>(syst)->BaseName());
60  }
61 
62  // Open output file
63  TFile* f = TFile::Open("covmxmanager.root", "read");
64 
65  // Get and save matrix
66  CovMxManager* mgr = CovMxManager::LoadFrom(f).release();
67  covmx::CovarianceMatrix* mx = mgr->GetCovarianceMatrix(samples, systs);
68  mx->SaveTo(TFile::Open("textmx.root", "recreate"));
69 
70  // Draw full matrix
71  TH2* h = mx->GetFullCovMxTH2();
72  TCanvas c("c", "c", 1600, 900);
73  h->Draw("colz");
74  c.SaveAs("plots/covmxmgr/fullmx.png");
75 
76  // Draw oscillated matrix
77  auto calc = DefaultSterileCalc(4);
79  mx->Predict(samples, calc);
80  std::vector<double> pred;
81  for (covmx::Sample& sample : samples) {
82  TH1D* hTmp = sample.GetPrediction()->Predict(calc).ToTH1(sample.GetPOT());
83  for (size_t i = 1; i <= (size_t)hTmp->GetNbinsX(); ++i) {
84  pred.push_back(hTmp->GetBinContent(i));
85  }
86  }
87  h = mx->GetCovMxRelativeTH2(pred);
88  h->Draw("colz");
89  c.SaveAs("plots/covmxmgr/oscmx.png");
90 
91  h->GetZaxis()->SetRangeUser(-0.1, 0.3);
92  c.SaveAs("plots/covmxmgr/oscmxzoom.png");
93 
94  h = mx->GetCovMxAbsoluteTH2();
95  h->Draw("colz");
96  c.SetLogz(true);
97  h->GetZaxis()->SetRangeUser(0.01, 1e9);
98  c.SaveAs("plots/covmxmgr/oscmxabs.png");
99  c.SetLogz(false);
100 
101  h = mx->GetCorrelationMxTH2();
102  h->Draw("colz");
103  c.SaveAs("plots/covmxmgr/corrmx.png");
104 
105  for (std::string syst : systs) {
106  mx = mgr->GetCovarianceMatrix(samples, {syst});
107  mx->Predict(samples, calc);
108  h = mx->GetCovMxRelativeTH2(pred);
109  h->Draw("colz");
110  if (h->GetMinimum() == 0) h->SetMinimum(-1e-3*h->GetMaximum());
111  c.SaveAs(Form("plots/covmxmgr/indiv/%s.png", syst.c_str()));
112  }
113 
114  } catch(const std::exception& e) {
115  std::cout << std::endl << "Exception occurred! " << e.what() << std::endl;
116  }
117 
118 } // macro MakeCovMx
119 
TH2D * GetCovMxRelativeTH2(std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SetPrediction(covmx::Sample &sample, bool systs=true)
Definition: Utilities.h:564
const std::vector< std::string > fdComps
Definition: TestCovMxNew.C:18
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
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
#define pot
const std::vector< std::string > ndComps
Definition: TestCovMxNew.C:12
std::vector< std::string > comps
OStream cout
Definition: OStream.cxx:6
std::vector< const ISyst * > getAllXsecSysts_2018_RPAFix()
TDirectory * dir
Definition: macro.C:5
static std::unique_ptr< CovMxManager > LoadFrom(TDirectory *dir)
Standard implementation.
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
std::vector< const ISyst * > LoadSysts(covmx::Sample sample)
Get systematics for a given sample.
Definition: Utilities.h:524
std::vector< TH1D * > GetSpectraFromDir(TDirectory *dir, std::vector< std::string > comps, double pot)
virtual void SaveTo(TDirectory *dir, const std::string &name) const
const double kAna2018FHCPOT
Definition: Exposures.h:207
Eigen::MatrixXd Predict(std::vector< Sample > samples, osc::IOscCalcAdjustable *calc, const SystShifts &systs=SystShifts::Nominal())
Float_t e
Definition: plot.C:35
void TestCovMxNew(std::string opt)
Definition: TestCovMxNew.C:34
const double kAna2018FHCLivetime
Definition: Exposures.h:213
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
enum BeamMode string