GenieInterface.cxx
Go to the documentation of this file.
1 /*
2  * GenieInterface.cxx:
3  * Interface to GENIE event records.
4  *
5  * Created on: Oct. 19, 2018
6  * Author: J. Wolcott <jwolcott@fnal.gov>
7  */
8 
9 #include "NOvARwgt/interfaces/GenieInterface.h"
10 
11 #include <memory>
12 #include <regex>
13 #include <unordered_map>
14 
15 #include "TParticle.h"
16 
17 
18 // GENIE includes
19 #if GENIE_MAJOR_VERSION >= 3
23 #else
24 #include "Conventions/Constants.h"
25 #include "Conventions/GVersion.h"
26 #include "EVGCore/EventRecord.h"
27 #include "GHEP/GHepParticle.h"
28 #include "ReWeight/GSystUncertainty.h"
29 #include "ReWeight/GReWeight.h"
30 #include "ReWeight/GReWeightNuXSecCCQE.h"
31 #endif
32 
33 #include "NOvARwgt/rwgt/EventRecord.h"
34 #include "NOvARwgt/util/GeneratorSupportConfig.h"
35 
36 // not in the public interface, so not in NOvARwgt/inc
37 #include "../rwgt/genie/GenieInternalTools.h"
38 
39 namespace novarwgt
40 {
41  // todo: need to decide how to specify which weights should be calculated on-the-fly if they're not pre-stored
42  novarwgt::EventRecord ConvertGenieEvent(const genie::EventRecord *evt,
43  const ReweightList &storedWgts)
44  {
45  novarwgt::EventRecord rec;
46 
47  // assume that the version of GENIE used to read the event was also the one used to make it.
48  // not foolproof, but best we can do here.
49  rec.generator = kGENIE;
50  rec.generatorVersion = novarwgt::internal::GetGENIEVersion();
51  rec.generatorConfigStr = novarwgt::internal::GetGENIETune();
52 
53  // BEWARE: there are two ways to calculate many of the values below:
54  // (1) using the 'selected' kinematics,
55  // i.e., those as seen by the off-shell hit nucleon inside the nucleus;
56  // (2) the kinematics as you'd calculate using the final-state particles.
57  // They can be significantly different.
58  // The values stored inside NOvASoft
59  // (and therefore the ones used to construct many of the weights)
60  // correspond to scenario (2).
61  //
62  // The calculations here are cribbed from nutools/GENIE/GENIEHelper.cxx,
63  // which is the 'authoritative' source of the values we have
64  // (it's how they enter the event stream in NOvASoft).
65  // Please see commentary there on the calculations
66  bool hitnucl = evt->Summary()->InitState().Tgt().HitNucIsSet();
67  TLorentzVector & k1 = *((evt->Probe())->P4());
68  TLorentzVector & k2 = *((evt->FinalStatePrimaryLepton())->P4());
69  TLorentzVector q = k1-k2; // q=k1-k2, 4-p transfer
70 
71  double Q2 = -1 * q.M2(); // momemtum transfer
72  double v = q.Energy(); // v (E transfer to the had system)
73  double y = v/k1.Energy(); // Inelasticity, y = q*P1/k1*P1
74  double W2, W;
75  W2 = W = -1;
76  if ( hitnucl || evt->Summary()->ProcInfo().IsCoherent() )
77  {
78  const double M = genie::constants::kNucleonMass;
79 // x = 0.5*Q2/(M*v);
80  W2 = M*M + 2*M*v - Q2;
81  W = std::sqrt(W2);
82  }
83 
84  rec.nupdg = evt->Probe()->Pdg();
85  rec.isCC = evt->Summary()->ProcInfo().IsWeakCC();
86  rec.reaction = novarwgt::ReactionType(evt->Summary()->ProcInfo().ScatteringTypeId()); // these are the same enum so just cast
87  rec.struckNucl = evt->Summary()->InitState().Tgt().HitNucPdg();
88 
89  rec.Enu = evt->Probe()->E();
90  rec.y = y;
91  rec.W = W;
92  rec.q = q;
93 
94  rec.A = evt->Summary()->InitState().Tgt().A();
95 
96  rec.npiplus = evt->NEntries(211, genie::kIStHadronInTheNucleus);
97  rec.npizero = evt->NEntries(111, genie::kIStHadronInTheNucleus);
98  rec.npiminus = evt->NEntries(-211, genie::kIStHadronInTheNucleus);
99 
100  for (const TObject * obj : (*evt))
101  {
102  auto particle = dynamic_cast<const genie::GHepParticle*>(obj);
103  if (particle->Status() != genie::kIStStableFinalState)
104  continue;
105 
106  rec.fsPartMult[particle->Pdg()]++; // default-constructor initializes to 0
107  rec.fsPartKE[particle->Pdg()] += particle->KinE(); // ditto
108  }
109 
110  rec.genieWeights = storedWgts;
111 
112  rec.origGenieEvt = evt;
113 
114  return rec;
115  }
116 
117 }
bool IsWeakCC(void) const
double E(void) const
Get energy.
Definition: GHepParticle.h:92
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
int HitNucPdg(void) const
Definition: Target.cxx:321
int A(void) const
Definition: Target.h:71
T sqrt(T number)
Definition: d0nt_math.hpp:156
caf::StandardRecord * rec
Definition: tutCAFMacro.C:20
virtual unsigned int NEntries(int pdg, GHepStatus_t ist, int start=0) const
Definition: GHepRecord.cxx:513
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
int Pdg(void) const
Definition: GHepParticle.h:64
int evt
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:370
novarwgt::EventRecord ConvertGenieEvent(const genie::EventRecord *evt, const ReweightList &storedWgts)
ScatteringType_t ScatteringTypeId(void) const
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
bool HitNucIsSet(void) const
Definition: Target.cxx:300
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
const Target & Tgt(void) const
Definition: InitialState.h:67
#define W(x)
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:40
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97