novarwgt_sr_interface_test.cxx
Go to the documentation of this file.
1 /*
2  * novarwgt_sr_interface_test.cxx
3  *
4  * Test the NOvA StandardRecord interface to NOvARwgt.
5  *
6  * Created on: Jan. 18, 2018
7  * Author: J. Wolcott <jwolcott@fnal.gov>
8  */
9 
10 #include <memory>
11 #include <unordered_map>
12 
13 #include "TDatabasePDG.h"
14 
15 #include "NOvARwgt/rwgt/EventRecord.h"
16 #include "NOvARwgt/test/TestEventTools.h"
17 
21 
22 
23 namespace novarwgt
24 {
25  namespace test
26  {
27  const std::unordered_map<novarwgt::ReactionType, caf::mode_type_> REACTION_ENUM_MAP
28  {
42  };
43 
44  const std::unordered_map<novarwgt::Generator, caf::generator_> GENERATOR_ID_MAP
45  {
46  {novarwgt::kGENIE, caf::kGENIE}
47  };
48 
49  /// Utility function to convert an EventRecord back into a StandardRecord.
50  caf::StandardRecord ToStandardRecord(const novarwgt::EventRecord &rec)
51  {
52  caf::StandardRecord srrec;
53  caf::SRNeutrino nu;
54 
55  nu.generator = GENERATOR_ID_MAP.find(rec.generator) != GENERATOR_ID_MAP.end()
56  ? GENERATOR_ID_MAP.at(rec.generator)
58  nu.genVersion.insert(nu.genVersion.begin(), rec.generatorVersion.begin(), rec.generatorVersion.end());
59  nu.genConfigString = rec.generatorConfigStr;
60 
61  nu.pdg = rec.nupdg;
62  nu.iscc = rec.isCC;
63  nu.mode = REACTION_ENUM_MAP.at(rec.reaction);
64  nu.hitnuc = rec.struckNucl;
65 
66  nu.E = rec.Enu;
67  nu.q2 = -rec.q * rec.q;
68  nu.y = rec.y;
69  nu.W2 = rec.W * rec.W;
70 
71  nu.tgtA = rec.A;
72 
73  nu.npiplus = rec.npiplus;
74  nu.npizero = rec.npizero;
75  nu.npiminus = rec.npiminus;
76 
77  // the final-state hadrons
78  for (const auto & fsMultPair : rec.fsPartMult)
79  {
80  // we'll divide it up equally amongst the particles since we don't save any more info
81  int pdg = fsMultPair.first;
82  const auto itKE = rec.fsPartKE.find(pdg);
83  double totalKE = ((itKE == rec.fsPartKE.end()) ? 0 : itKE->second);
84  unsigned int nPart = fsMultPair.second;
85  double KEPerParticle = totalKE / nPart;
86  double mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
87  double EPerParticle = KEPerParticle + mass;
88  double MomPerParticle = sqrt(EPerParticle*EPerParticle - mass*mass);
89  for (unsigned int i = 0; i < nPart; i++)
90  {
92  part.pdg = pdg;
93  part.p = caf::SRLorentzVector();
94  part.p.E = EPerParticle;
95  part.p.px = MomPerParticle;
96  part.p.py = 0;
97  part.p.pz = 0;
98  nu.prim.emplace_back(std::move(part));
99  }
100  }
101 
102 
103  nu.rwgt.genie.resize(rec.genieWeights.size());
104  for (std::size_t idx = 0; idx < rec.genieWeights.size(); idx++)
105  {
106  if (rec.genieWeights.IsSet(idx))
107  {
108  auto nugenIdx = std::find_if(novarwgt::kNugenKnobTranslationTable.begin(),
110  [idx](const auto & pair)
111  {
112  return std::size_t(pair.second) == idx;
113  })->first;
114 
115  nu.rwgt.genie[nugenIdx].minus2sigma = rec.genieWeights[idx].minus2sigma;
116  nu.rwgt.genie[nugenIdx].minus1sigma = rec.genieWeights[idx].minus1sigma;
117  nu.rwgt.genie[nugenIdx].plus1sigma = rec.genieWeights[idx].plus1sigma;
118  nu.rwgt.genie[nugenIdx].plus2sigma = rec.genieWeights[idx].plus2sigma;
119  }
120  }
121 
122  nu.isvtxcont = !rec.expectNoWeights;
123 
124  srrec.mc.nu.emplace_back(std::move(nu));
125  srrec.mc.nnu = 1;
126 
127  return std::move(srrec);
128  }
129  }
130 }
131 
132 
133 int main()
134 {
135  std::cout << "NOvARwgt self-test: NOvA StandardRecord (caf::StandardRecord)" << std::endl;
136  std::cout << "=============================================================" << std::endl;
137  std::cout << std::endl;
138 
139  novarwgt::InputVals params
140  {
141  {"EmpiricalMEC", true},
142  };
143 
144  auto testEventsIn = novarwgt::test::GetTestEvents();
145  decltype(testEventsIn) testEventsOut;
146  for (const auto & evtInPair : testEventsIn)
147  {
148  std::cout << "Converting event: '" << evtInPair.first << "'" << std::endl;
149 
150  const auto & evt = evtInPair.second;
151  // this may look a little silly
152  // (converting novarwgt::EventRecords to the StandardRecord types and then back)
153  // but (assuming the ToStandardRecord() function above is written correctly)
154  // it tests that ConvertSREvent() is working correctly...
155  auto sr = novarwgt::test::ToStandardRecord(evt.Event());
156  auto newEvt = novarwgt::ConvertSREvent(&sr, true);
157 
158  // since I can't copy the derived class pointer without knowing its type, I have to dynamic_cast it first. ugh
159  std::unique_ptr<std::exception> exc = nullptr;
160  if (evt.ExpectedException())
161  exc = std::make_unique<novarwgt::UnsupportedGeneratorException>(*dynamic_cast<const novarwgt::UnsupportedGeneratorException*>(evt.ExpectedException()));
162 
163  testEventsOut.emplace( evtInPair.first,
164  novarwgt::test::TestEvent<novarwgt::EventRecord>(newEvt, evt.ExpectedWeight(), evt.Tune(), evt.ExpectedSysts(), std::move(exc)));
165  }
166  bool ok = novarwgt::test::Validate(testEventsOut, params, true);
167  if (ok)
168  std::cout << "All " << testEventsOut.size() << " events produced the expected behavior." << std::endl;
169 
170  return ok ? 0 : 1;
171 }
unsigned int npizero
Number of &#39;s after neutrino reaction, before FSI.
Definition: SRNeutrino.h:65
float q2
Squared momentum transfer [GeV^2].
Definition: SRNeutrino.h:77
caf::StandardRecord ToStandardRecord(const novarwgt::EventRecord &rec)
Utility function to convert an EventRecord back into a StandardRecord.
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< SRTrueParticle > prim
Primary daughters, lepton comes first in vector.
Definition: SRNeutrino.h:84
unsigned int npiplus
Number of &#39;s after neutrino reaction, before FSI.
Definition: SRNeutrino.h:63
const novarwgt::EventRecord & ConvertSREvent(const T *sr, bool forceNoCache)
Copy information out of a StandardRecord.
SRMCReweight rwgt
Definition: SRNeutrino.h:91
std::vector< unsigned int > genVersion
Version of the generator that created this neutrino interaction.
Definition: SRNeutrino.h:88
caf::StandardRecord * rec
Definition: tutCAFMacro.C:20
const std::unordered_map< rwgt::ReweightLabel_t, novarwgt::ReweightKnob > kNugenKnobTranslationTable
int tgtA
A of target nucleus.
Definition: SRNeutrino.h:75
std::string genConfigString
String associated with generator configuration. (For GENIE 3, this is the "Comprehensive Model Config...
Definition: SRNeutrino.h:89
The SRNeutrino is a representation of neutrino interaction information.
Definition: SRNeutrino.h:21
int main()
MAIN FUNCTION.
unsigned int npiminus
Number of &#39;s after neutrino reaction, before FSI.
Definition: SRNeutrino.h:64
TString part[npart]
Definition: Style.C:32
bool iscc
true if charged-current interaction, false if not.
Definition: SRNeutrino.h:60
float W2
Invariant mass of final state squared. [GeV^2].
Definition: SRNeutrino.h:80
const std::unordered_map< novarwgt::ReactionType, caf::mode_type_ > REACTION_ENUM_MAP
int evt
caf::StandardRecord * sr
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
int mode
interaction mode from enum mode_type::[QE, RES, COH, ...]
Definition: SRNeutrino.h:59
4-vector with more efficient storage than TLorentzVector
OStream cout
Definition: OStream.cxx:6
short pdg
pdg code
Definition: SRNeutrino.h:25
bool isvtxcont
Checks if neutrino true vertex is within detector.
Definition: SRNeutrino.h:38
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
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
float E
True energy [GeV].
Definition: SRNeutrino.h:28
std::vector< SRGenieWeights > genie
GENIE weights.
Definition: SRMCReweight.h:23
float y
Bjorken y = (p.q) / (k.p), fractional energy loss of incoming particle [Dimensionless].
Definition: SRNeutrino.h:79
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
SRLorentzVector p
Momentum 4-vector.
generator_ generator
The generator that created this neutrino interaction.
Definition: SRNeutrino.h:87