EHadVisMECpairs.C
Go to the documentation of this file.
1 /*
2  * EHadVisMECpairs.C:
3  * What's the EhadVis response of different nucleon-nucleon pair types in Empirical MEC?
4  *
5  * Created on: Apr 13, 2018
6  * Author: J. Wolcott <jwolcott@fnal.gov>
7  */
8 
9 #include "TCanvas.h"
10 #include "TH1.h"
11 #include "TLatex.h"
12 #include "TLegend.h"
13 
15 #include "CAFAna/Core/Loaders.h"
16 #include "CAFAna/Core/Spectrum.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
18 #include "CAFAna/Cuts/TruthCuts.h"
21 #include "CAFAna/Vars/TruthVars.h"
22 #include "CAFAna/Vars/XsecTunes.h"
24 
25 namespace jw
26 {
27  const std::map<int, std::string> NN_PAIR_NAMES
28  {
29  { 2000000200, "nn" },
30  { 2000000201, "np" },
31  { 2000000202, "pp" },
32  };
33 }
34 
35 void EHadVisMECpairs(bool makeSpecs=true, std::string fileName="", bool makePlots=true)
36 {
37  std::map<ana::Loaders::FluxType, std::map<int, std::unique_ptr<ana::Spectrum>>> nnPairSpecs;
38 
39  if (makeSpecs)
40  {
41  std::map<ana::Loaders::FluxType, ana::Prod4NomLoaders> loadersMap;
42  for (const auto & hc : {ana::Loaders::kRHC, ana::Loaders::kFHC})
43  loadersMap.emplace(std::piecewise_construct,
44  std::forward_as_tuple(hc),
45  std::forward_as_tuple(ana::kNumuConcat, hc));
46  for (auto & loadersPair : loadersMap)
47  loadersPair.second.SetSpillCut(ana::kStandardSpillCuts);
48 
49  ana::HistAxis ehadVisAxis {"Visible E_{had} (GeV)", ana::Binning::Simple(75, 0, 0.5), ana::kNumuHadVisE};
50 
51  for (auto & loadersPair : loadersMap)
52  {
53  ana::Loaders::FluxType fluxType = loadersPair.first;
55  auto & loaders = loadersPair.second;
56  for (const auto & NNpair : jw::NN_PAIR_NAMES)
57  {
58  if ( (fluxType == ana::Loaders::kRHC && NNpair.second == "nn")
59  || (fluxType == ana::Loaders::kFHC && NNpair.second == "pp") )
60  continue;
61 
62  nnPairSpecs[fluxType].emplace(NNpair.first,
63  std::make_unique<ana::Spectrum>(
64  loaders.GetLoader(caf::kNEARDET, ana::Loaders::kMC),
65  ehadVisAxis,
66  nuCut && ana::kIsDytmanMEC && (ana::kHitNuc == NNpair.first),
69  )
70  );
71  }
72  }
73 
74  for (auto & loadersPair : loadersMap)
75  loadersPair.second.Go();
76 
77  if (!fileName.empty())
78  {
79 
80  TFile outf (fileName.c_str (), "recreate");
81  if (outf.IsZombie ())
82  {
83  std::cerr << "Couldn't open output file: " << fileName << std::endl;
84  std::cerr << "Abort." << std::endl;
85  return;
86  }
87 
88  for (auto & specsPair : nnPairSpecs)
89  {
90  ana::Loaders::FluxType fluxType = specsPair.first;
91 
92  auto fluxDir = outf.mkdir(fluxType == ana::Loaders::kFHC ? "FHC" : "RHC");
93  for (const auto & specPair : specsPair.second)
94  {
95  TDirectory * specDir = fluxDir->mkdir(jw::NN_PAIR_NAMES.at(specPair.first).c_str());
96  specPair.second->SaveTo (specDir);
97  }
98  } // for (specsPair)
99  } // if (!fileName.empty())
100  } // if (makeSpecs)
101  else if (!fileName.empty())
102  {
103  TFile inf(fileName.c_str());
104  if (inf.IsZombie ())
105  {
106  std::cerr << "Couldn't open input file: " << fileName << std::endl;
107  std::cerr << "Abort." << std::endl;
108  return;
109  }
110  for (const auto & fluxType : {ana::Loaders::kFHC, ana::Loaders::kRHC})
111  {
112  auto * fluxDir = dynamic_cast<TDirectory*>(inf.Get(fluxType == ana::Loaders::kFHC ? "FHC" : "RHC"));
113  if (!fluxDir)
114  continue;
115 
116  for (const auto & nnPairPair : jw::NN_PAIR_NAMES)
117  {
118  TDirectory * dir = dynamic_cast<TDirectory*>(fluxDir->Get(nnPairPair.second.c_str()));
119  if (!dir)
120  continue;
121 
122  nnPairSpecs[fluxType][nnPairPair.first] = ana::Spectrum::LoadFrom(dir);
123 // std::cout << " loaded spectrum for " << (fluxType == ana::Loaders::kFHC ? "FHC" : "RHC")
124 // << " " << nnPairPair.second << std::endl;
125  }
126  } // for (specsPair)
127 
128 
129  }
130  else
131  {
132  std::cerr << "Not making or loading spectra. Nothing to do!" << std::endl;
133  return;
134  }
135 
136  if (!makePlots)
137  {
138  std::cout << "Not requested to make plots. All done." << std::endl;
139  return;
140  }
141 
142  std::cout << "Making plots..." << std::endl;
143 
144  for (const auto & fluxPair : nnPairSpecs)
145  {
146  TCanvas c("", "", 600, 500);
147  c.SetLeftMargin(0.15);
148  c.SetBottomMargin(0.15);
149  TLegend leg(0.6, 0.7, 0.9, 0.9);
150  leg.SetFillStyle(0);
151  leg.SetBorderSize(0);
152  TH1 * hfirst = nullptr;
153  double ymax = 0;
154  for (const auto & specPair : fluxPair.second)
155  {
156  int nnPairType = specPair.first;
157  const ana::Spectrum & spec = *(specPair.second);
158  TH1 * h = spec.ToTH1(spec.POT(), nnPairType == 2000000201 ? kRed : kBlue);
159  h->Scale(1./h->Integral());
160  h->GetXaxis()->SetTitle("Visible E_{had} (GeV)");
161  h->GetXaxis()->CenterTitle();
162  h->GetXaxis()->SetTitleOffset(1.25);
163  h->GetYaxis()->SetTitle("Fraction of events");
164  h->GetYaxis()->CenterTitle();
165  h->GetYaxis()->SetTitleOffset(1.25);
166  if (!hfirst)
167  {
168  h->Draw("][ hist");
169  hfirst = h;
170  }
171  else
172  h->Draw("][ hist same");
173  leg.AddEntry(h, (jw::NN_PAIR_NAMES.at(nnPairType) + " initial state").c_str(), "l");
174  ymax = std::max(ymax, h->GetMaximum());
175  }
176  hfirst->GetYaxis()->SetRangeUser(0, ymax*1.2);
177  leg.Draw();
178  std::string fluxName = fluxPair.first == ana::Loaders::kRHC ? "nubar" : "nu";
179  TLatex fluxLabel(0.2, 0.83, fluxPair.first == ana::Loaders::kRHC ? "Antineutrino beam" : "Neutrino beam");
180  fluxLabel.SetNDC();
181  fluxLabel.SetTextSize(0.05);
182  fluxLabel.Draw();
183  util::SaveObj(&c, "EhadVis_MEC_nnPairs_" + fluxName,
184  "/nova/ana/users/jwolcott/2p2h_tuning/2018/plots",
185  {".png", ".pdf", ".root"});
186  }
187 }
Near Detector underground.
Definition: SREnums.h:10
T max(const caf::Proxy< T > &a, T b)
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
fileName
Definition: plotROC.py:78
enum BeamMode kRed
TFile * inf
Definition: drawXsec.C:1
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:148
void makePlots()
Definition: makePlots.C:17
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
OStream cerr
Definition: OStream.cxx:7
const std::map< int, std::string > NN_PAIR_NAMES
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
const Var kHitNuc
Definition: TruthVars.h:19
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
Double_t ymax
Definition: plot.C:25
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
TFile * outf
Definition: testXsec.C:51
double POT() const
Definition: Spectrum.h:227
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
TDirectory * dir
Definition: macro.C:5
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void EHadVisMECpairs(bool makeSpecs=true, std::string fileName="", bool makePlots=true)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
void SaveObj(const TObject *obj, const std::string &filenameStub, const std::string &dirName="", const std::vector< std::string > exts={".png",".eps",".root"}, bool silent=false)
Definition: ROOTHelpers.h:21
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
enum BeamMode string