MakePPFXRatios.C
Go to the documentation of this file.
1 
2 #include "TString.h"
3 #include "TFile.h"
4 
6 
7 #include "CAFAna/Core/ISyst.h"
8 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Core/Sample.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
11 #include "CAFAna/Vars/HistAxes.h"
12 #include "CAFAna/Vars/XsecTunes.h"
15 
16 #include "NuXAna/Cuts/NusCuts18.h"
18 
19 using namespace ana;
20 
21 void MakePPFXRatios(TString opt) {
22 
23  covmx::Sample sample = opt.Contains("neardet", TString::kIgnoreCase) ?
25 
26  // Here get loaders
28  if (sample.detector == covmx::kNearDet) loaders.SetND(true);
30 
31  // This should surely be a vector of loaders and cuts, right?
32  auto det = (sample.detector == covmx::kNearDet) ? caf::kNEARDET : caf::kFARDET;
33 
34  std::vector<std::tuple<SpectrumLoaderBase*, const Cut, std::string>> components;
35 
36  SpectrumLoaderBase* nonswap = &loaders.GetLoader(det,
38 
39  // Unoscillated components (for both detectors)
40  components.push_back(std::make_tuple(nonswap, kIsBeamNue && !kIsAntiNu, "cc_nue_to_nue" ));
41  components.push_back(std::make_tuple(nonswap, kIsBeamNue && kIsAntiNu, "cc_nuebar_to_nuebar" ));
42  components.push_back(std::make_tuple(nonswap, kIsNumuCC && !kIsAntiNu, "cc_numu_to_numu" ));
43  components.push_back(std::make_tuple(nonswap, kIsNumuCC && kIsAntiNu, "cc_numubar_to_numubar"));
44  components.push_back(std::make_tuple(nonswap, kIsNC, "nc" ));
45 
46  // CC appearance channels (for far detector only)
47  if (sample.detector == covmx::kFarDet) {
50  SpectrumLoaderBase* tau = &loaders.GetLoader(det,
52  components.push_back(std::make_tuple(swap, kIsSig && !kIsAntiNu, "cc_numu_to_nue" ));
53  components.push_back(std::make_tuple(swap, kIsSig && kIsAntiNu, "cc_numubar_to_nuebar" ));
54  components.push_back(std::make_tuple(swap, kIsNumuApp && !kIsAntiNu, "cc_nue_to_numu" ));
55  components.push_back(std::make_tuple(swap, kIsNumuApp && kIsAntiNu, "cc_nuebar_to_numubar" ));
56  components.push_back(std::make_tuple(tau, kIsTauFromE && !kIsAntiNu, "cc_nue_to_nutau" ));
57  components.push_back(std::make_tuple(tau, kIsTauFromE && kIsAntiNu, "cc_nuebar_to_nutaubar" ));
58  components.push_back(std::make_tuple(tau, kIsTauFromMu && !kIsAntiNu, "cc_numu_to_nutau" ));
59  components.push_back(std::make_tuple(tau, kIsTauFromMu && kIsAntiNu, "cc_numubar_to_nutaubar"));
60  }
61 
62  // PPFX weight vars
63  std::vector<Var> kPPFXFluxUnivWgts;
64  for (size_t i = 0; i < 100; ++i) {
65  const Var tempPPFXWgt([i](const caf::SRProxy* sr){
66  if (!sr->hdr.ismc) return 1.f;
67  if (sr->mc.nnu != 1) return 1.f;
68  if (sr->mc.nu[0].rwgt.ppfx.vuniv[i] <= 0) return 1.f;
69  return (float)sr->mc.nu[0].rwgt.ppfx.vuniv[i];
70  });
71  kPPFXFluxUnivWgts.push_back(tempPPFXWgt);
72  }
73 
74  // Analysis cuts and vars
75  const HistAxis axis = sample.detector == covmx::kNearDet ? kNCNDAxisE : kNCFDAxisE;
76  const Cut sel = sample.detector == covmx::kNearDet ? kNus18ND : kNus18FD;
77  const Var kReweight = kXSecCVWgt2018RPAFix * kPPFXFluxCVWgt;
78 
79  // Initialise maps
80  std::vector<Spectrum*> nominal(components.size());
81  std::vector<std::vector<Spectrum*> > ppfx(components.size(),
82  std::vector<Spectrum*>(kPPFXFluxUnivWgts.size(), nullptr));
83 
84  // Loop over each component
85  for (size_t c = 0; c < components.size(); ++c) {
86 
87  // Nominal spectra
88  nominal[c] = new Spectrum(*std::get<0>(components[c]),
89  axis, sel && std::get<1>(components[c]), kNoShift, kReweight);
90 
91  // PPFX universe spectra
92  for (size_t u = 0; u < 100; ++u)
93  ppfx[c][u] = new Spectrum(*std::get<0>(components[c]),
94  axis, sel && std::get<1>(components[c]), kNoShift, kXSecCVWgt2018RPAFix * kPPFXFluxUnivWgts[u]);
95 
96  } // for component
97 
98  loaders.Go();
99 
100  // Save em all
101  std::string fileName = "ppfx_ratios_" + sample.GetTag() + ".root";
102  TFile* f = TFile::Open(fileName.c_str(), "recreate");
103  TDirectory* outDir = f->mkdir(sample.GetTag().c_str());
104 
105  TDirectory* dir = outDir->mkdir("nominal");
106  for (size_t c = 0; c < components.size(); ++c) {
107  nominal[c]->SaveTo(dir->mkdir(std::get<2>(components[c]).c_str()));
108  }
109 
110  dir = outDir->mkdir("ppfx");
111  for (size_t u = 0; u < 100; ++u) {
112  TDirectory* udir = dir->mkdir(Form("u%zu",u+1));
113  for (size_t c = 0; c < components.size(); ++c) {
114  ppfx[c][u]->SaveTo(udir->mkdir(std::get<2>(components[c]).c_str()));
115  }
116  }
117 
118  delete f;
119 
120 }
121 
Near Detector underground.
Definition: SREnums.h:10
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
fileName
Definition: plotROC.py:78
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
std::vector< double > Spectrum
Definition: Constants.h:746
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
const Cut kIsBeamNue(CCFlavSel(12, 12))
Select CC .
const Cut kIsNumuApp(CCFlavSel(14, 12))
Select CC .
std::string GetTag() const
Definition: Sample.cxx:94
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const Cut kNus18FD
Definition: NusCuts18.h:100
std::string outDir
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
const Cut kIsTauFromE(CCFlavSel(16, 12))
Select CC .
const Cut kIsTauFromMu(CCFlavSel(16, 14))
Select CC .
if(dump)
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
Detector detector
Definition: Sample.h:97
const Cut kIsSig(CCFlavSel(12, 14))
Select CC .
caf::StandardRecord * sr
void MakePPFXRatios(std::string detector, int job, bool concat=true, bool rhc=false)
const Var kXSecCVWgt2018RPAFix
Definition: XsecTunes.h:52
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
const SystShifts kNoShift
Definition: SystShifts.cxx:22
Base class for the various types of spectrum loader.
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
TDirectory * dir
Definition: macro.C:5
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void SetND(bool nd)
Definition: Loaders.h:65
caf::Proxy< bool > ismc
Definition: SRProxy.h:242
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
enum BeamMode string