nutools_test.cxx
Go to the documentation of this file.
1 /*
2  * nutools_test.cxx
3  *
4  * Test the NuTools interface.
5  *
6  * Created on: Jan. 17, 2018
7  * Author: J. Wolcott <jwolcott@fnal.gov>
8  */
9 
10 #include <iostream>
11 #include <memory>
12 
13 #include "TDatabasePDG.h"
14 
15 #include "NOvARwgt/rwgt/EventRecord.h"
16 #include "NOvARwgt/interfaces/NuToolsInterface.h"
17 #include "NOvARwgt/test/TestEventTools.h"
18 
23 
24 
25 namespace novarwgt
26 {
27  namespace test
28  {
29  struct NuSimDataEvt
30  {
31  std::unique_ptr<simb::MCParticle> nuPart;
32  std::unique_ptr<simb::MCParticle> lepPart;
33  std::vector<std::unique_ptr<simb::MCParticle>> hadParts;
34 
35  std::unique_ptr<simb::MCTruth> mctruth;
36  std::unique_ptr<simb::GTruth> gtruth;
37  };
38 
39  /// Utility function to convert an EventRecord back into NuSimData components.
40  NuSimDataEvt ToNuSimData(const novarwgt::EventRecord &rec)
41  {
43 
44  evt.mctruth = std::make_unique<simb::MCTruth>();
45 
46  std::unordered_map<std::string, std::string> config;
47  if (!rec.generatorConfigStr.empty())
48  config["tune"] = rec.generatorConfigStr;
49  evt.mctruth->SetGeneratorInfo(rec.generator == novarwgt::kGENIE ? simb::Generator_t::kGENIE : simb::Generator_t::kUnknown,
50  EncodeGeneratorVersion(rec.generatorVersion),
51  config);
52 
53  TLorentzVector fourMom;
54  fourMom.SetE(rec.Enu);
56  evt.nuPart = std::make_unique<simb::MCParticle>(1, rec.nupdg, "");
57  evt.nuPart->AddTrajectoryPoint({0,0,0,0}, fourMom);
58  evt.mctruth->Add(*evt.nuPart);
59 
60  fourMom.SetE((1-rec.y)*rec.Enu);
61  evt.lepPart = std::make_unique<simb::MCParticle>(2, rec.isCC ? (rec.nupdg > 0 ? rec.nupdg-1 : rec.nupdg+1) : rec.nupdg, "", 1);
62  evt.lepPart->AddTrajectoryPoint({0,0,0,0}, fourMom);
63  evt.mctruth->Add(*evt.lepPart);
64 
65  // not exactly right, but we aren't currently storing Z
66  // format is +10LZZZAAAI
67  int targetPdg = (rec.A == 1) ? 1000010010 : 1000000000 + (rec.A/2)*10000 + rec.A*10;
68  evt.mctruth->SetNeutrino(
69  rec.isCC ? simb::kCC : simb::kNC,
70  simb::int_type_(rec.reaction),
71  simb::int_type_(rec.reaction+1000), // not exactly right, but we don't use this
72  targetPdg,
73  rec.struckNucl,
74  0, // not storing struck quark
75  rec.W,
76  0, // not storing x
77  rec.y,
78  rec.Q2());
79 
80  evt.gtruth = std::make_unique<simb::GTruth>();
81  evt.gtruth->fNumPi0 = rec.npizero;
82  evt.gtruth->fNumPiMinus = rec.npiminus;
83  evt.gtruth->fNumPiPlus = rec.npiplus;
84 
85  // the final-state hadrons
86  for (const auto & fsMultPair : rec.fsPartMult)
87  {
88  // we'll divide it up equally amongst the particles since we don't save any more info
89  int pdg = fsMultPair.first;
90  const auto itKE = rec.fsPartKE.find(pdg);
91  double totalKE = ((itKE == rec.fsPartKE.end()) ? 0 : itKE->second);
92  unsigned int nPart = fsMultPair.second;
93  double KEPerParticle = totalKE / nPart;
94  double mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
95  double EPerParticle = KEPerParticle + mass;
96  double MomPerParticle = sqrt(EPerParticle*EPerParticle - mass*mass);
97  for (unsigned int i = 0; i < nPart; i++)
98  {
99  int particleIdx = 2+i+1; // 2=(original neutrino + lepton above)
100  evt.hadParts.emplace_back(std::make_unique<simb::MCParticle>(particleIdx, pdg, "", 1, mass));
101  TLorentzVector fourMom(MomPerParticle, 0, 0, EPerParticle);
102  evt.hadParts.back()->AddTrajectoryPoint({0,0,0,0}, fourMom);
103  evt.mctruth->Add(*evt.hadParts.back());
104 
105  }
106  }
107 
108 
109  return evt;
110  }
111  }
112 }
113 
114 
115 int main()
116 {
117  std::cout << "NOvARwgt self-test: NuTools (simb::MCTruth and simb::GTruth)" << std::endl;
118  std::cout << "============================================================" << std::endl;
119  std::cout << std::endl;
120 
121  novarwgt::InputVals params
122  {
123  {"EmpiricalMEC", true},
124  };
125 
126  auto testEventsIn = novarwgt::test::GetTestEvents();
127  std::cout << "Considering " << testEventsIn.size() << " events:" << std::endl;
128  decltype(testEventsIn) testEventsOut;
129  for (const auto & evtInPair : testEventsIn)
130  {
131  const auto & evt = evtInPair.second;
132  // this may look a little silly
133  // (converting novarwgt::EventRecords to the NuSimData types and then back)
134  // but (assuming the ToNuSimData() function above is written correctly)
135  // it tests the ConvertNuToolsEvent() is working correctly...
136  auto nuSimDataEvt = novarwgt::test::ToNuSimData(evt.Event());
137  auto newEvt = novarwgt::ConvertNuToolsEvent(nuSimDataEvt.mctruth.get(), nuSimDataEvt.gtruth.get(), evt.Event().genieWeights);
138 
139  // since I can't copy the derived class pointer without knowing its type, I have to dynamic_cast it first. ugh
140  std::unique_ptr<std::exception> exc = nullptr;
141  if (evt.ExpectedException())
142  exc = std::make_unique<novarwgt::UnsupportedGeneratorException>(*dynamic_cast<const novarwgt::UnsupportedGeneratorException*>(evt.ExpectedException()));
143 
144  testEventsOut.emplace( evtInPair.first,
145  novarwgt::test::TestEvent<novarwgt::EventRecord>(newEvt,
146  evt.ExpectedWeight(),
147  evt.Tune(),
148  evt.ExpectedSysts(),
149  std::move(exc)));
150  }
151  bool ok = novarwgt::test::Validate(testEventsOut, params, true);
152  if (ok)
153  std::cout << "All " << testEventsOut.size() << " events produced the expected behavior." << std::endl;
154 
155  return ok ? 0 : 1;
156 }
std::unique_ptr< simb::MCTruth > mctruth
std::unique_ptr< simb::MCParticle > lepPart
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::StandardRecord * rec
Definition: tutCAFMacro.C:20
std::unique_ptr< simb::MCParticle > nuPart
Particle class.
std::unique_ptr< simb::GTruth > gtruth
novarwgt::EventRecord ConvertNuToolsEvent(const simb::MCTruth *mctruth, const simb::GTruth *gtruth, const novarwgt::ReweightList &rwList)
NuSimDataEvt ToNuSimData(const novarwgt::EventRecord &rec)
Utility function to convert an EventRecord back into NuSimData components.
int evt
std::vector< std::unique_ptr< simb::MCParticle > > hadParts
std::string EncodeGeneratorVersion(const std::vector< int > &verVec)
OStream cout
Definition: OStream.cxx:6
int main()
int_type_
Neutrino interaction categories.
Definition: MCNeutrino.h:79