12 #include "CAFAna/Core/Binning.h" 18 #include "CAFAna/Cuts/NumuCuts.h" 23 #include "CAFAna/Vars/NumuVars.h" 45 TLatex*
prelim =
new TLatex(.9, .95,
"NO#nuA Preliminary");
46 prelim->SetTextColor(
kBlue);
48 prelim->SetTextSize(2/30.);
49 prelim->SetTextAlign(32);
56 if(
sr->
mc.
nnu == 0)
return false;
114 std::cout <<
"\nrun : ---- Running CAF NuMU ND Validation.\n";
123 TH1D* spillHist =
new TH1D(
"reco-spills",
";Detector;Spills", 3, 0, 3);
141 std::vector<Selection> selections_base;
142 std::vector<Selection> selections_int;
143 std::vector<Selection> selections_contain;
146 selections_base.emplace_back(
kIsNumuCC,
"true_numucc",
147 "Selected events are truly numu cc");
148 if(isFD) selections_base.emplace_back(
kNumuFD,
"numuFD",
149 "Selected events pass the numu FD selection");
150 else selections_base.emplace_back(
kNumuND,
"numuND",
151 "Selected events pass the numu ND selection");
153 selections_int.emplace_back(
kNoCut,
"no_cut",
154 "No selection applied");
155 selections_int.emplace_back(
kIsQE,
"true_qe",
156 "Selected events are truly quasi-elastic");
157 selections_int.emplace_back(!
kIsQE,
"true_non_qe",
158 "Selected events are truly non quasi-elastic");
160 selections_contain.emplace_back(
kNoCut,
"no_cut",
161 "No selection applied");
164 "Selected events are truly contained in the ND");
166 "Selected events are truly un-contained in the ND");
172 const Var kTrueNHitSlc({
"mc.nu",
"mc.nnu",
"mc.nu.nhitslc"},
175 const Var kTrueNHitTot({
"mc.nu",
"mc.nnu",
"mc.nu.nhittot"},
178 const Var kRecoNHitTot({
"slc.nhit"},
181 const Var kTruehadE({
"mc.nu",
"mc.nnu",
"mc.nu.*"},
183 {
return (
sr->
mc.
nnu == 0 ||
sr->
mc.
nu[0].prim.size() == 0) ? 0 :
186 const Var kTruehadESlc({
"mc.nu",
"mc.nnu",
"mc.nu.*"},
188 {
return (
sr->
mc.
nnu == 0 ||
sr->
mc.
nu[0].prim.size() == 0) ? 0 :
189 (
sr->
mc.
nu[0].visEinslc -
sr->
mc.
nu[0].prim[0].visEinslc);
191 const Var kRecohadE({
"energy.numusimp.hadcalE"},
193 {
return (
sr ->
energy.numusimp.hadcalE);
197 {
return (
sr->
mc.
nu[0].prim.size() == 0) ? 0 :
198 (
sr->
mc.
nu[0].prim[0].p.E());
200 const Var kTrueMuESlc({
"mc.nu",
"mc.nnu",
"mc.nu.*"},
202 {
return (
sr->
mc.
nu[0].prim.size() == 0) ? 0 :
203 sr->
mc.
nu[0].prim[0].p.E();
205 const Var kRecoMuE({
"energy.mutrkE.E",
"sel.remid.bestidx"},
214 {
return (
sr->
mc.
nnu == 0) ? 0 :
217 const Var kRecoNuEnergy({
"energy.numusimp.trkccE"},
219 {
return (
sr->
mc.
nnu == 0) ? 0 :
225 const Var kRMTOT_NHitSlc({
"mc.nu",
"mc.nnu",
"mc.nu.nhitslc",
"slc.nhit"},
228 if(
sr->
mc.
nnu == 0)
return -10.;
231 double nhitTrue = (
double)
sr->
mc.
nu[0].nhitslc;
232 return ((nhitReco - nhitTrue) / nhitTrue);
237 const Var kRMTOT_Energy({
"mc.nu",
"mc.nnu",
"mc.nu.E",
"energy.numusimp.trkccE"},
239 {
return (
sr->
mc.
nnu == 0) ? -10 :
243 const Var kRMTOT_hadE({
"mc.nu",
"mc.nnu",
"mc.nu.*",
"energy.numusimp.hadcalE"},
246 if(
sr->
mc.
nnu == 0)
return -10.f;
247 if(
sr->
mc.
nu[0].prim.size() == 0)
return -10.
f;
248 float trueHadE = (
sr ->
mc.nu[0].E) - (
sr ->
mc.nu[0].prim[0].p.E());
250 float recoHadE=
sr->
energy.numusimp.trkccE-muE;
252 if(trueHadE == 0)
return -10.f;
253 return ( (recoHadE - trueHadE) / trueHadE );
256 const Var kRMTOT_Emu({
"mc.nu",
"mc.nnu",
"mc.nu.*",
"energy.mutrkE.E",
"sel.remid.bestidx"},
259 if(
sr->
mc.
nnu == 0)
return -10.f;
260 if(
sr->
mc.
nu[0].prim.size() == 0)
return -10.
f;
262 if(
sr->
sel.
remid.bestidx == 999)
return -10.f;
263 float trueMuE = (
sr ->
mc.nu[0].prim[0].p.E());
265 if(trueMuE == 0)
return -10.f;
266 return ( (recoMuE - trueMuE) / trueMuE );
281 const HistAxis axRMTOT_Energy(
"Neutrino energy (Reco. - true)/true",
Binning::Simple(100,-2,2), kRMTOT_Energy);
288 const HistAxis axTrueHadE(
"Had E. (GeV)", kHadronicEnergyBinning, kTruehadE);
289 const HistAxis axTrueHadESlc(
"Had E. (GeV)", kHadronicEnergyBinning, kTruehadESlc);
290 const HistAxis axRecoHadE(
"Had E. (GeV)", kHadronicEnergyBinning, kRecohadE);
292 const HistAxis axTrueMuESlc(
"Muon E. (GeV)", kEnergyBinning, kTrueMuESlc);
295 const HistAxis axRecoNuEnergy(
"Neutrino E. (GeV)", kEnergyBinning, kRecoNuEnergy);
299 std::vector<TangibleAxis> variables;
300 variables.emplace_back(axRMTOT_NHitSlc,
"res-nhit_slc_reco_minus_true_over_true",
301 "Number of hits in slc. (reco-true)/true");
302 variables.emplace_back(axRMTOT_Energy,
"res-neutrino_energy_reco_minus_true_over_true",
303 "Neutrino energy (reco-true)/true");
304 variables.emplace_back(axRMTOT_HadE,
"res-hadronic_energy_reco_minus_true_over_true",
305 "Hadronic energy (reco-true)/true");
306 variables.emplace_back(axRMTOT_Emu,
"res-muon_energy_reco_minus_true_over_true",
307 "Muon energy (reco-true)/true");
309 std::vector<TangibleAxis> dist_variables;
325 std::vector<UsefulHist>
hists;
326 hists.reserve(selections_base.size() * selections_int.size() * selections_contain.size() * variables.size() + selections_base.size() * selections_int.size() * selections_contain.size() * dist_variables.size());
330 for(
const auto& sel_base:selections_base){
331 for(
const auto& sel_int:selections_int){
332 for(
const auto& sel_contain:selections_contain){
333 for(
const auto& variable:variables){
334 hists.emplace_back(sel_base, sel_int, sel_contain, variable, loader);
346 for(
const auto& sel_base:selections_base){
347 for(
const auto& sel_int:selections_int){
348 for(
const auto& variable:variables){
349 hists.emplace_back(sel_base, sel_int, selections_contain[0], variable, loader);
358 std::cout <<
"\nrun : --- run loader.\n";
361 double pot = hists[0].fHist.POT();
363 std::cout <<
"\nrun : --- save output.\n";
365 TFile outputfile(kOutputFileName.c_str(),
"RECREATE");
367 for(
const auto&
hist:hists){
370 temp->Write(tempName.c_str());
375 TH1F *pothist =
new TH1F(
"meta-TotalPOT",
"TotalPOT", 1, 0, 1);
376 TH1F *evthist =
new TH1F(
"meta-TotalEvents",
"TotalEvents", 1, 0, 1);
377 pothist->SetBinContent(1, pot);
378 evthist->SetBinContent(1, spillHist->GetEntries());
379 pothist->Write(
"meta-TotalPOT");
380 evthist->Write(
"meta-TotalEvents");
382 std::cout <<
"\nrun : --- close output files.\n";
const Cut kTruthContainedND([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(fabs(sr->mc.nu[0].vtx.X()) > 200|| fabs(sr->mc.nu[0].vtx.Y()) > 200|| sr->mc.nu[0].vtx.Z()< 0|| sr->mc.nu[0].vtx.Z() > 1650) return false;return true;})
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
Tangible< Cut > Selection
void caf_numu_reco_minus_true(std::string kInputFileName, std::string kOutputFileName, bool isFD)
void SetSpillCut(const SpillCut &cut)
Representation of a spectrum in any variable, with associated POT.
virtual void AddSpillHistogram(TH1 *h, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
Uses include counting the total POT or spills in a run.
Tangible(const T &obj, const std::string &shortName, const std::string &blurb)
const Binning kHadronicEnergyBinning
SRRemid remid
Output from RecoMuonID (ReMId) package.
Tangible< HistAxis > TangibleAxis
const Binning kEnergyBinning
virtual void Go() override
Load all the registered spectra.
short nnu
Number of neutrinos in nu vector (0 or 1)
unsigned int nhit
number of hits
The StandardRecord is the primary top-level object in the Common Analysis File trees.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
UsefulHist(Selection selA, Selection selB, Selection selC, TangibleAxis tanAxis, SpectrumLoader &loader, const Cut &bkg=kNoCut)
SRIDBranch sel
Selector (PID) branch.
assert(nhit_max >=nhit_nbins)
void Preliminary()
Put NOvA Preliminary on plots.
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
SRSlice slc
Slice branch: nhit, extents, time, etc.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
Template for Var and SpillVar.
SREnergyBranch energy
Energy estimator branch.
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.