NOvARwgtInterface.cxx
Go to the documentation of this file.
1 /*
2  * NOvARwgtInterface.cxx:
3  * Create NOvARwgt objects from StandardRecords
4  *
5  * Created on: Jan. 24, 2019
6  * Author: J. Wolcott <jwolcott@fnal.gov>
7  */
8 
10 
11 #include <unordered_map>
12 #include <unordered_set>
13 
14 #include "nugen/NuReweight/ReweightLabels.h"
15 
19 
20 #include "NOvARwgt/rwgt/EventRecord.h"
21 
22 // use novarwgt namespace because this is effectively a NOvARwgt file
23 // (just ported here so that whenever we update SR this doesn't get out of date)
24 namespace novarwgt
25 {
26  // declared here, not in the .h, so that it doesn't escape this file
27  namespace
28  {
29  const std::unordered_map<rwgt::EReweightLabel, novarwgt::ReweightKnob> RW_KNOB_NAME_MAP
30  {
31  {rwgt::fReweightMaCCQE, novarwgt::kKnob_MaCCQE},
32  {rwgt::fReweightMaCCRES, novarwgt::kKnob_MaCCRES},
33  };
34 
35  // there's some slight skew between the ordering of these in the GENIE & SR versions
36  const std::unordered_map<caf::mode_type_, novarwgt::ReactionType> REACTION_ENUM_MAP
37  {
46 
51 
53 
54  // these two don't have GENIE versions...
55  //{caf::kEM, novarwgt::kScNull},
56  //{caf::kWeakMix, novarwgt::kScNull },
57 
58  };
59 
60  const std::unordered_map<caf::generator_, novarwgt::Generator> GENERATOR_ID_MAP
61  {
62  {caf::kGENIE, novarwgt::kGENIE}
63  };
64 
65  // End-of-program warn-er
66  class UnmatchedKnobWarn
67  {
68  public:
69  UnmatchedKnobWarn() = default;
70 
71  ~UnmatchedKnobWarn()
72  {
73  if (fUnmatchedKnobs.empty())
74  return;
75 
76  std::cerr << "WARNING: no corresponding NOvARwgt knob found for the following rwgt::ReweightLabel indices (found in events):" << std::endl;
77  for (const auto & knob : fUnmatchedKnobs)
78  std::cerr << " " << knob << ",";
80  std::cerr << "Stored weights for them were ignored." << std::endl;
81  }
82 
83  void AddKnob(std::size_t newKnob)
84  {
85  fUnmatchedKnobs.insert(newKnob);
86  }
87 
88  private:
89  std::unordered_set<std::size_t> fUnmatchedKnobs;
90  };
91  }
92 
93  /// Copy information out of an SRNeutrino
94  template <typename T>
95  const novarwgt::EventRecord& ConvertSRTruth(const T * nu, bool forceNoCache)
96  {
98  "novarwgt::ConvertSRTruth() can only be used with SRNeutrino or SRNeutrinoProxy");
99 
100  // hold a cache of the last StandardRecord we converted to a novarwgt::EventRecord
101  // so that syst operations that use the same event repeatedly
102  // don't convert it lots of times in a row
103  static novarwgt::EventRecord sCachedEvt;
104 
105  bool matches = std::abs(nu->E - sCachedEvt.Enu) < 0.0001 && std::abs(nu->q2 - sCachedEvt.Q2()) < 0.0001
106  && std::abs(nu->y - sCachedEvt.y) < 0.0001 && (nu->iscc == sCachedEvt.isCC);
107 
108  if (matches && !forceNoCache)
109  return sCachedEvt;
110 
111  sCachedEvt.Reset();
112 
113  sCachedEvt.generator = GENERATOR_ID_MAP.at(nu->generator);
114  for (const auto & v : nu->genVersion)
115  sCachedEvt.generatorVersion.push_back(v);
116  sCachedEvt.generatorConfigStr = nu->genConfigString;
117 
118  sCachedEvt.nupdg = nu->pdg;
119  sCachedEvt.isCC = nu->iscc;
120  sCachedEvt.reaction = REACTION_ENUM_MAP.at(caf::mode_type_(int(nu->mode))); // nu->mode was originally a caf::mode_type_
121  sCachedEvt.struckNucl = nu->hitnuc;
122 
123  sCachedEvt.Enu = nu->E;
124  sCachedEvt.q = TLorentzVector(1, 0, 0, nu->y * nu->E); // need to make up the direction since we didn't store it
125  sCachedEvt.q.SetRho(sqrt(nu->q2 + sCachedEvt.q.E()*sCachedEvt.q.E()));
126  sCachedEvt.y = nu->y;
127  sCachedEvt.W = sqrt(nu->W2);
128 
129  sCachedEvt.A = nu->tgtA;
130 
131  sCachedEvt.npiplus = nu->npiplus;
132  sCachedEvt.npizero = nu->npizero;
133  sCachedEvt.npiminus = nu->npiminus;
134 
135  for (const auto & part : nu->prim)
136  {
137  sCachedEvt.fsPartMult[part.pdg]++; // default constructor initializes to 0
138  sCachedEvt.fsPartKE[part.pdg] += (1 - 1./part.p.Gamma()) * part.p.E;
139  }
140 
141  // the following doesn't work because SRMCReweightProxy owns a VectorProxy<SRGenieWeightsProxy>,
142  // and a VectorProxy fundamentally contains a vector<T*> instead of a vector<T>.
143 // static_assert(sizeof(novarwgt::ReweightVals) == sizeof(caf::SRGenieWeights),
144 // "The NOvARwgt-StandardRecord interface only works if novarwgt::ReweightVals and caf::SRGenieWeights have the same size!");
145 // std::vector<novarwgt::ReweightVals> systWgts(nu->rwgt.genie.size());
146 // std::memcpy(systWgts.data(), nu->rwgt.genie.data(), nu->rwgt.genie.size() * sizeof(novarwgt::ReweightVals));
147 // sCachedEvt.genieWeights = ReweightList(std::move(systWgts));
148 
149  static UnmatchedKnobWarn warn;
150  // since we have to dereference the elements of the vector to copy them anyway,
151  // we may as well just copy them one by one.
152  sCachedEvt.genieWeights.resize(novarwgt::kLastKnob);
153  for (std::size_t idx = 0; idx < nu->rwgt.genie.size(); idx++)
154  {
155  // slightly less typing
156  const auto & stdRecWgts = nu->rwgt.genie[idx];
157  auto rwgtIdx = static_cast<rwgt::ReweightLabel_t>(idx);
158  try
159  {
160  auto novaRwIdx = novarwgt::kNugenKnobTranslationTable.at(rwgtIdx);
161  sCachedEvt.genieWeights[novaRwIdx]
162  = {stdRecWgts.minus2sigma, stdRecWgts.minus1sigma, stdRecWgts.plus1sigma, stdRecWgts.plus2sigma};
163  }
164  catch(std::out_of_range & e)
165  {
166  warn.AddKnob(rwgtIdx);
167  }
168  }
169 
170  if (!nu->isvtxcont)
171  sCachedEvt.expectNoWeights = true;
172 
173  // GENIE reweight table maker doesn't make weights for nu-on-e events with Q^2 < 0
174  // since the GENIE reweight calcs tend to flip out
175  if (nu->rwgt.genie.empty() && nu->mode == caf::kElectronScattering)
176  sCachedEvt.expectNoWeights = true;
177 
178  return sCachedEvt;
179 
180  }
181  template const novarwgt::EventRecord& ConvertSRTruth(const caf::SRNeutrino *, bool);
182  template const novarwgt::EventRecord& ConvertSRTruth(const caf::SRNeutrinoProxy *, bool);
183 
184  /// Copy information out of a StandardRecord
185  template <typename T>
186  const novarwgt::EventRecord& ConvertSREvent(const T * sr, bool forceNoCache)
187  {
189  "novarwgt::ConvertSRTruth() can only be used with SRNeutrino or SRNeutrinoProxy");
190 
191  static novarwgt::EventRecord emptyRecord;
192  if (sr->mc.nnu != 1)
193  return emptyRecord;
194 
195  return ConvertSRTruth( &(sr->mc.nu[0]), forceNoCache );
196 
197 
198  }
199  template const novarwgt::EventRecord& ConvertSREvent(const caf::StandardRecord *, bool);
200  template const novarwgt::EventRecord& ConvertSREvent(const caf::SRProxy *, bool);
201 }
mode_type_
Neutrino interaction categories.
Definition: SREnums.h:47
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
T sqrt(T number)
Definition: d0nt_math.hpp:156
const novarwgt::EventRecord & ConvertSREvent(const T *sr, bool forceNoCache)
Copy information out of a StandardRecord.
OStream cerr
Definition: OStream.cxx:7
const std::unordered_map< rwgt::ReweightLabel_t, novarwgt::ReweightKnob > kNugenKnobTranslationTable
float abs(float number)
Definition: d0nt_math.hpp:39
The SRNeutrino is a representation of neutrino interaction information.
Definition: SRNeutrino.h:21
std::unordered_set< std::size_t > fUnmatchedKnobs
TString part[npart]
Definition: Style.C:32
const XML_Char int const XML_Char * value
Definition: expat.h:331
const std::unordered_map< novarwgt::ReactionType, caf::mode_type_ > REACTION_ENUM_MAP
caf::StandardRecord * sr
def warn(msg)
print a warning message, from: https://cdcvs.fnal.gov/redmine/projects/novaart/repository/entry/trunk...
Definition: common_tools.py:16
void out_of_range(const char *function, int max, int index, const char *msg1="", const char *msg2="")
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const std::unordered_map< novarwgt::Generator, caf::generator_ > GENERATOR_ID_MAP
double T
Definition: Xdiff_gwt.C:5
const novarwgt::EventRecord & ConvertSRTruth(const T *nu, bool forceNoCache)
Copy information out of an SRNeutrino.
Float_t e
Definition: plot.C:35