KinematicsHistos.C
Go to the documentation of this file.
1 /*
2  * KinematicsHistos.C
3  *
4  * Created on: Feb. 15, 2018
5  * (copied from 2017 version)
6  * Author: J. Wolcott <jwolcott@fnal.gov>
7  */
8 
9 
10 #ifdef __CINT__
11 #include <iostream>
12 void KinematicsHistos(std::string hornCurrent, std::string productionName, std::string dataMC = "")
13 {
14  std::cerr << "CINT not supported (or tolerated) here. Please compile." << std::endl;
15  return 1;
16 }
17 
18 #else
19 
20 #include <algorithm>
21 #include <exception>
22 #include <iomanip>
23 #include <memory>
24 #include <unordered_map>
25 #include <set>
26 #include <string>
27 
28 #include "TFile.h"
29 #include "TMath.h"
30 #include "TParameter.h"
31 
34 
35 #include "CAFAna/Analysis/SALoaders.h"
37 #include "CAFAna/Cuts/Cuts.h"
39 #include "CAFAna/Cuts/SpillCuts.h"
40 #include "CAFAna/Cuts/TruthCuts.h"
41 #include "CAFAna/Core/Spectrum.h"
43 #include "CAFAna/Vars/Vars.h"
47 #include "CAFAna/Vars/TruthVars.h"
49 
50 
51 // /////////////////////////////////////
52 //
53 // custom types for use below
54 
55 const std::set<std::string> KNOWN_DATASETS
56 {
57  "test",
58  "prod3",
59  "prod4",
60 };
61 
62 const std::string BASE_DIR = "/nova/ana/users/jwolcott/2p2h_tuning/2018/hists";
63 
64 // I don't want to 'use' the whole namespace, but I want this function.
65 // (allows joining of requirements with '+' operator.)
66 using ana::operator+;
67 
68 ////////
69 
70 
71 const ana::Var kVisMuHadE
72  ([](const caf::SRProxy * sr)
73  {
74  return kMuE(sr) + ana::kNumuHadVisE(sr);
75  }
76 );
77 
78 
79 ////////
80 
81 struct IntConfig
82 {
83  IntConfig(const std::string& nm, const ana::Cut & cut, unsigned int color)
84  : name(nm), selCut(cut), histColor(color) {};
87  unsigned int histColor;
88 };
89 
90 ana::Cut Q3CutFactory(float lb, float ub)
91 {
92  return std::move(ana::Cut([lb, ub](const caf::SRProxy * sr)
93  {
94  double q3 = ana::kRecoQmag(sr);
95  return lb <= q3 && q3 <= ub;
96  }
97  ));
98 }
99 // /////////////////////////////////////
100 //
101 // configuration
102 const std::unordered_map<std::string, ana::HistAxis> OBSERVABLE_HIST_AXES {
103  {"Ehad", {"Reco E_{had} (GeV)", ana::Binning::Simple(50, 0., 2.), kHadE}},
104  {"EhadVis", {"Visible E_{had} (GeV)", ana::Binning::Simple(50, 0., 1.), ana::kNumuHadVisE}},
105  {"Emu", {"Reco E_{#mu} (GeV)", ana::Binning::Simple(50, 0., 3.), kMuE}},
106  {"Enu", {"Reco E_{#nu} (GeV)", ana::Binning::Simple(50, 0., 4.), kCCE}},
107  {"Emu+EHadVis", {"Reco E_{#mu} + Vis. E_{had} (GeV)", ana::Binning::Simple(50, 0., 4.), kVisMuHadE}},
108  {"RecoQ0", {"Reco 'q_{0}' = E_{had,vis} (GeV)", ana::Binning::Simple(50, 0., 2.), ana::kNumuHadVisE}},
109  {"RecoQ2", {"Reco Q^{2} (GeV^{2})", ana::Binning::Simple(50, 0., 2.), ana::kRecoQ2}},
110  {"RecoQmag", {"Reco |#vec{q}| = #sqrt{Q^{2}+q_{0}^{2}} (GeV)", ana::Binning::Simple(50, 0., 2.), ana::kRecoQmag}},
111  {"RecoW", {"Reco W (GeV)", ana::Binning::Simple(50, 0., 2.), ana::kRecoW}},
112 };
113 
114 const std::unordered_map<std::string, ana::HistAxis> TRUTH_HIST_AXES {
115  {"TrueEnu", {"True E_{#nu} (GeV)", ana::Binning::Simple(50, 0., 4.), ana::kTrueQ0}},
116  {"TrueQ2", { "True Q^{2} (GeV^{2})", ana::Binning::Simple(50, 0., 2.), ana::kTrueQ2}},
117  {"TrueQmag", {"True |#vec{q}| (GeV)", ana::Binning::Simple(50, 0., 2.), ana::kTrueQ3}},
118  {"TrueQ0", {"True q_{0} (GeV)", ana::Binning::Simple(50, 0., 2.), ana::kTrueQ0}},
119  {"TrueW", {"True W (GeV)", ana::Binning::Simple(50, 0., 2.), ana::kTrueW}},
120 };
121 
122 
123 const std::unordered_map<caf::mode_type_, IntConfig, std::hash<int>> INT_CONFIG {
124  {caf::kQE, IntConfig("QE", ana::kIsNumuCC && ana::kIsQE, kAzure+7)},
126  {caf::kDIS, IntConfig("DIS", ana::kIsNumuCC && ana::kIsDIS, kGray)},
128 
130  IntConfig("Other", !ana::kIsNumuCC || !(ana::kIsQE || ana::kIsRes || ana::kIsDIS || ana::kIsDytmanMEC), kBlack) }
131 };
132 
137 
138 //kRescaleMAQE_NT --ok2
139 //* kRPAWeightCCQE2018_ST --ok2
140 //* kRPAWeightRES2017_ST --ok1
141 //* kFixNonres1PiST --ok1
142 //* kRescaleHighWDIS_NT --ok1
143 //* kEmpiricalMECWgt2018_NT --ok2
144 
145 const std::unordered_map<std::string, ana::Var> MECQE_WGTS {
148  { "MINERvA-1p1h", ana::kMINERvA_Wgt_QE}, // note that MINERvA weight INCLUDES RPA
149  { "MINERvA-2p2h", ana::kMINERvA_Wgt_MEC},
150  { "MINERvA-2p2hnp", ana::kMINERvA_Wgt_MECNP},
151  { "MINERvA-2p2hpp", ana::kMINERvA_Wgt_MECPP},
155 };
156 
157 // ugh.
158 const std::unordered_map<std::string, std::unordered_map<std::string, ana::Cut>> SLICE_BREAKDOWNS {
159  { "RecoQ0", // will replace the plot in the list above with sliced. so don't call it 'Ehad' or 'EhadVis'
160  {
161  {"+q3=0-0.1", Q3CutFactory(0, 0.1) },
162  {"+q3=0.1-0.2", Q3CutFactory(0.1, 0.2) },
163  {"+q3=0.2-0.3", Q3CutFactory(0.2, 0.3) },
164  {"+q3=0.3-0.4", Q3CutFactory(0.3, 0.4) },
165  {"+q3=0.4-0.5", Q3CutFactory(0.4, 0.5) },
166  {"+q3=0.5-0.6", Q3CutFactory(0.5, 0.6) },
167  {"+q3=0.6-0.7", Q3CutFactory(0.6, 0.7) },
168  {"+q3=0.7-0.8", Q3CutFactory(0.7, 0.8) },
169  {"+q3=0.8-0.9", Q3CutFactory(0.8, 0.9) },
170  {"+q3=0.9-1.0", Q3CutFactory(0.9, 1.0) },
171  }
172  }
173 
174 };
175 
176 // /////////////////////////////////////
177 
179 {
180  public:
181  IOManager(std::string hornCurrent, std::string dataSet, std::string dataMC)
182  : fDataSet(dataSet), fHornCurrent(hornCurrent), fDataMC(ana::Loaders::kMC),
183  fLoader(nullptr)
184  {
185  if (dataMC.size() > 0)
186  {
187  std::transform(dataMC.begin(), dataMC.end(), dataMC.begin(), ::tolower);
188  if (dataMC == "mc")
189  fDataMC = ana::Loaders::kMC;
190  else if (dataMC == "data")
191  fDataMC = ana::Loaders::kData;
192  else
193  throw std::runtime_error("Choose 'data' or 'mc'. You passed: '" + dataMC + "'");
194 
195  }
196 
197  if (hornCurrent.size() > 0)
198  {
199  std::transform(hornCurrent.begin(), hornCurrent.end(), hornCurrent.begin(), ::toupper);
200  if (hornCurrent != "FHC" && hornCurrent != "RHC")
201  throw std::runtime_error("Choose 'FHC' or 'RHC'. You passed: '" + hornCurrent + "'");
202  fHornCurrent = hornCurrent;
203  }
204  }
205 
206  ana::Loaders::DataMC GetDataMC() const { return fDataMC; };
207  ana::SpectrumLoaderBase & GetLoader() const;
208  std::string TagNameBySample(const std::string & name) const;
209 
210  private:
214 
215  mutable std::unique_ptr<ana::SpectrumLoaderBase> fLoader;
216 
217 };
218 
220 {
221  // cached after the first time
222  if (fLoader != nullptr)
223  return *fLoader;
224 
225 
226  if (KNOWN_DATASETS.find(fDataSet) == KNOWN_DATASETS.end())
227  throw std::runtime_error("Unrecognized dataset: '" + fDataSet + "'");
228 
229 
231  if (fHornCurrent == "FHC")
232  {
233  if (fDataSet == "prod4")
234  switch (fDataMC)
235  {
236  case ana::Loaders::kData:
237  fname = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns_numu2018";
238  break;
239 
240  case ana::Loaders::kMC:
241  fname = "prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018";
242  break;
243  }
244  else if (fDataSet == "test")
245  {
246  switch(fDataMC)
247  {
248  case ana::Loaders::kMC:
249  fname = "/pnfs/nova/persistent/production/concat/R17-11-14-prod4reco.d/nd/sumdecaf/numu_contain/genie/*period5*57*.root";
250  break;
251 
252  case ana::Loaders::kData:
253  fname = "/pnfs/nova/persistent/production/concat/R17-09-05-prod4recopreview.f/nd/sumdecaf/numu_contain/numi/*period5*7*.root";
254  break;
255  }
256  } // if (dataSet == "test")
257  } // if (FHC)
258  else if (fHornCurrent == "RHC")
259  {
260  if (fDataSet == "prod4")
261  switch (fDataMC)
262  {
263  case ana::Loaders::kData:
264  fname = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns_numu2018";
265  break;
266 
267  case ana::Loaders::kMC:
268  fname = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018";
269  break;
270  }
271  else if (fDataSet == "test")
272  {
273  switch(fDataMC)
274  {
275  case ana::Loaders::kMC:
276  fname = "/pnfs/nova/persistent/production/concat/R17-11-14-prod4reco.e/nd/sumdecaf/numu_contain/genie/*period4*57*.root";
277  break;
278 
279  case ana::Loaders::kData:
280  fname = "/pnfs/nova/persistent/production/concat/R17-09-05-prod4recopreview.f/nd/sumdecaf/numu_contain/numi/*period4*7*.root";
281  break;
282  }
283  } // if (dataSet == "test")
284  }
285 
286  if (!fLoader && fname != "")
287  fLoader.reset(new ana::SpectrumLoader(fname));
288 
289  fLoader->SetSpillCut(ana::kStandardSpillCuts);
290  return *fLoader;
291 }
292 
293 ////////
294 
296 {
297  // split at last period, if any
298  std::string work(name);
299  if (name.find("/") != std::string::npos)
300  work = name.substr(name.rfind("/"));
301  auto splitPoint = work.rfind('.');
302  std::string basename(work.substr(0, splitPoint));
303  std::string ext(name.substr(splitPoint));
304 
305  std::string infix = fDataSet + "_" + fHornCurrent + "_";
306  switch (fDataMC)
307  {
308  case ana::Loaders::kMC:
309  infix += "MC";
310  break;
311 
312  case ana::Loaders::kData:
313  infix += "Data";
314  break;
315  }
316 
317  // on the grid we don't look for an absolute directory
318  std::string base_dir = getenv("_CONDOR_SCRATCH_DIR") ? "" : BASE_DIR + "/";
319  return std::move(base_dir + basename + "_" + infix + ext);
320 }
321 
322 // /////////////////////////////////////
323 //
324 // helper functions
325 
326 
327 ////////
328 
329 /// Plots of 'observables' (or things near to it): Emu, Ehad, Enu
330 /// If data, one Spectrum for each Var
331 /// If MC, one Spectrum per Var per interaction class, split further by which weights were applied
332 void MakeObservablesSpectra(const IOManager & mgr, std::unordered_map<std::string, std::unique_ptr<ana::Spectrum>> & spectra)
333 {
335  decltype(SLICE_BREAKDOWNS)::mapped_type no_slice { {"", ana::kNoCut} };
336 
337  for (const auto & axisPair : OBSERVABLE_HIST_AXES)
338  {
339  const auto & slices = (SLICE_BREAKDOWNS.find(axisPair.first) == SLICE_BREAKDOWNS.end()) ? no_slice : SLICE_BREAKDOWNS.at(axisPair.first);
340  for (const auto slicePair : slices)
341  {
342  auto plotname = axisPair.first + slicePair.first;
343  auto & axis = axisPair.second;
344  if (mgr.GetDataMC() == ana::Loaders::kData)
345  {
346  spectra.emplace(plotname, // var name
347  std::make_unique<ana::Spectrum>(mgr.GetLoader(), axis, selCut && slicePair.second)
348  );
349  continue;
350  }
351 
352  // MC is more complicated... breakdowns by interaction type; versions for various weightings
353  for (const auto & mecWgtPair : MECQE_WGTS)
354  {
355  auto wgts = mecWgtPair.second * kNonQEMECWgt;
356  for (const auto & intPair : INT_CONFIG)
357  {
358  std::string name = Form("%s{group=mecWgt,cat=%s}{group=interaction,cat=%s}", plotname.c_str(), mecWgtPair.first.c_str(), intPair.second.name.c_str());
359  // std::cout << " preparing Spectrum " << std::quoted(name) << std::endl;
360  auto cuts = selCut && slicePair.second && intPair.second.selCut;
361  spectra.emplace(name,
362  std::make_unique<ana::Spectrum>(mgr.GetLoader(), axis, cuts, ana::kNoShift, wgts)
363  );
364  } // for (intPair)
365  } // for (mecWgtPair)
366  } //for (slice)
367  } // for (configPair)
368 } // MakeObservableSpectra()
369 
370 ////////
371 
372 /// Plots of truth-only quantities (weights applied; true kinematics)
373 void MakeTruthSpectra(const IOManager & mgr, std::unordered_map<std::string, std::unique_ptr<ana::Spectrum>> & spectra)
374 {
375  for (const auto & intPair : INT_CONFIG)
376  {
377 
378  // plot the weights themselves
379  for (const auto & wgtPair : MECQE_WGTS)
380  {
381  std::string name = Form("weight{group=mecWgt,cat=%s}{group=interaction,cat=%s}", wgtPair.first.c_str(), intPair.second.name.c_str());
382  auto cuts = ana::kNumuND && intPair.second.selCut;
383  spectra.emplace(name,
384  std::make_unique<ana::Spectrum>(name, ana::Binning::Simple(100, 0, 2), mgr.GetLoader(), wgtPair.second, cuts)
385  );
386  } // for (wgtPair)
387 
388  // and true event spectra. useful for debugging what happened
389  for (const auto & mecWgtPair : MECQE_WGTS)
390  {
391  auto wgts = mecWgtPair.second * kNonQEMECWgt;
392  for (const auto & configPair : TRUTH_HIST_AXES)
393  {
394  const auto & name = configPair.first;
395  const auto & axis = configPair.second;
396  std::string plotName = Form("%s{group=mecWgt,cat=%s}{group=interaction,cat=%s}", configPair.first.c_str(), mecWgtPair.first.c_str(), intPair.second.name.c_str());
397  spectra.emplace(plotName,
398  std::make_unique<ana::Spectrum>(
399  mgr.GetLoader(),
400  axis,
401  intPair.second.selCut,
403  wgts
404  )
405  );
406 
407  } // for (config)
408 
409  } // for (xsecWgts)
410  } // for (intPair)
411 } // MakeTruthSpectra()
412 
413 
414 
415 // /////////////////////////////////////
416 
417 void SaveSpectra(std::unordered_map<std::string, std::unique_ptr<ana::Spectrum>> & spectra, const IOManager & mgr)
418 {
419  TFile outFile(mgr.TagNameBySample("rwgt_spectra.root").c_str(), "recreate");
420  std::cout << "Saving spectra to TFile: " << std::quoted(outFile.GetName()) << std::endl;
421  for (const auto & specPair : spectra)
422  specPair.second->SaveTo(&outFile, specPair.first.c_str());
423 }
424 
425 // /////////////////////////////////////
426 //
427 // main
428 
429 void KinematicsHistos(std::string hornCurrent, std::string productionName, std::string dataMC = "")
430 {
431  IOManager mgr(hornCurrent, productionName, dataMC);
432 
433  std::unordered_map<std::string, std::unique_ptr<ana::Spectrum>> allSpectra;
434  MakeObservablesSpectra(mgr, allSpectra);
435  if (mgr.GetDataMC() == ana::Loaders::kMC)
436  MakeTruthSpectra(mgr, allSpectra);
437 
438  // Do it!
439  mgr.GetLoader().Go();
440 
441  SaveSpectra(allSpectra, mgr);
442 
443 }
444 #endif // #if __CINT__
const Var kHadE
Definition: NumuVars.h:23
const Cut kIsQE
Definition: TruthCuts.cxx:104
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
const std::unordered_map< std::string, std::unordered_map< std::string, ana::Cut > > SLICE_BREAKDOWNS
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
void SaveSpectra(std::unordered_map< std::string, std::unique_ptr< ana::Spectrum >> &spectra, const IOManager &mgr)
const Var kMINERvA_Wgt_MEC
Definition: GenieWeights.h:232
const ana::Var kRecoQ2([](const caf::SRProxy *sr){const double M_mu_sqrd=util::sqr(0.1056);double E_mu=kMuE(sr);double p_mu=sqrt(util::sqr(E_mu)-M_mu_sqrd);return 2 *kCCE(sr)*(E_mu-p_mu *kCosNumi(sr))-M_mu_sqrd;})
Reconstructed four-momentum transfer invariant (Q^2)
Definition: NumuVars.h:146
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeObservablesSpectra(const IOManager &mgr, std::unordered_map< std::string, std::unique_ptr< ana::Spectrum >> &spectra)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
std::string fDataSet
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
Definition: NumuVars.h:149
const Var kEmpiricalMECtoValenciaMECWgt
Definition: GenieWeights.h:190
const Cut kIsRes
Definition: TruthCuts.cxx:111
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const std::unordered_map< caf::mode_type_, IntConfig, std::hash< int > > INT_CONFIG
const Color_t kMC
const std::set< std::string > KNOWN_DATASETS
IOManager(std::string hornCurrent, std::string dataSet, std::string dataMC)
static constexpr Double_t nm
Definition: Munits.h:133
std::string name
const Var kMINERvA_Wgt_MECPP
Definition: GenieWeights.h:240
const Var kRPAWeightRES2017
Definition: GenieWeights.h:115
OStream cerr
Definition: OStream.cxx:7
const ana::Var kNonQEMECWgt
ana::Cut selCut
const Cut kIsDIS
Definition: TruthCuts.cxx:118
const Var kTrueQ0
Definition: TruthVars.h:32
void KinematicsHistos(std::string hornCurrent, std::string productionName, std::string dataMC="")
const Var kMEC2018QElikeWgt
Definition: GenieWeights.h:159
const Var kEmpiricalMECWgt2018
See kEmpiricalMECWgt2018_NT.
Definition: GenieWeights.h:254
const std::unordered_map< std::string, ana::HistAxis > OBSERVABLE_HIST_AXES
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
ana::Cut Q3CutFactory(float lb, float ub)
const std::unordered_map< std::string, ana::Var > MECQE_WGTS
std::string fHornCurrent
const Var kTrueQ2
Definition: TruthVars.h:27
ana::Loaders::DataMC fDataMC
const Var kRPAWeightCCQE2018
Definition: GenieWeights.h:96
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
virtual void Go()=0
Load all the registered spectra.
void MakeTruthSpectra(const IOManager &mgr, std::unordered_map< std::string, std::unique_ptr< ana::Spectrum >> &spectra)
Plots of truth-only quantities (weights applied; true kinematics)
const Var kRescaleMAQE
Definition: GenieWeights.h:61
const Var kTrueQ3
Definition: TruthVars.h:38
TFile * outFile
Definition: PlotXSec.C:135
const ana::Var kVisMuHadE([](const caf::SRProxy *sr){return kMuE(sr)+ana::kNumuHadVisE(sr);})
unsigned int histColor
const Var kMINERvA_Wgt_QE
Definition: GenieWeights.h:228
const std::unordered_map< std::string, ana::HistAxis > TRUTH_HIST_AXES
std::string getenv(std::string const &name)
const Cut kNumuND
Definition: NumuCuts.h:55
const Var kCCE
Definition: NumuVars.h:21
caf::StandardRecord * sr
ana::Loaders::DataMC GetDataMC() const
IntConfig(const std::string &nm, const ana::Cut &cut, unsigned int color)
ana::SpectrumLoaderBase & GetLoader() const
const Var kMEC2018RESlikeWgt
Definition: GenieWeights.h:163
const Var kMINERvA_Wgt_MECNP
Definition: GenieWeights.h:236
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
Base class for the various types of spectrum loader.
static const double ub
Definition: Units.h:88
const Var kRecoW([](const caf::SRProxy *sr){const double M_p=0.938272;double WSq=M_p *M_p+2 *M_p *(kCCE(sr)-kMuE(sr))-kRecoQ2(sr);if(WSq > 0) return sqrt(WSq);else return-5.;})
Reconstructed invariant mass (W)
Definition: NumuVars.h:152
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
std::unique_ptr< ana::SpectrumLoaderBase > fLoader
std::string TagNameBySample(const std::string &name) const
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const Var kDytmanMEC_Disable
Definition: TruthVars.h:15
const std::string BASE_DIR
const Var kMuE
Definition: NumuVars.h:22
enum BeamMode kGreen
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kTrueW
Definition: TruthVars.h:22
const Var kRescaleHighWDIS
Definition: GenieWeights.h:316
const Var kFixNonres1Pi
Definition: GenieWeights.h:127
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
enum BeamMode string