Event library interface

Makefile

@@ -31,6 +31,7 @@ INITIAL_BUILD_TARGETS = print-make-info \
 		   physics-boosted-dark-matter \
 		   tools-flux-drivers \
 		   tools-geometry-drivers \
+		   tools-vmc \
 		   tools-masterclass
 FINAL_BUILD_TARGETS = doxygen-doc \
 		   apps \
@@ -202,6 +203,18 @@ else
 	@echo "** Building geometry-drivers was not enabled. Skipping..."
 endif
 
+tools-vmc: FORCE
+#ifeq ($(strip $(GOPT_ENABLE_VMC)),YES)
+	@echo " "
+	@echo "** Building VMC..."
+	cd ${GENIE}/src/Tools/VMC && \
+	$(MAKE) && \
+	cd ${GENIE}
+#else
+#	@echo " "
+#	@echo "** Building VMC was not enabled. Skipping..."
+#endif
+
 
 tools-masterclass: FORCE
 ifeq ($(strip $(GOPT_ENABLE_MASTERCLASS)),YES)
@@ -368,6 +381,7 @@ make-install-dirs: FORCE
 	mkdir ${GENIE_INC_INSTALLATION_PATH}/Tools/Flux
 	mkdir ${GENIE_INC_INSTALLATION_PATH}/Tools/Geometry
 	mkdir ${GENIE_INC_INSTALLATION_PATH}/Tools/Masterclass
+	mkdir ${GENIE_INC_INSTALLATION_PATH}/Tools/VMC
 
 
 copy-install-files: FORCE
@@ -424,6 +438,7 @@ copy-install-files: FORCE
 	cd ${GENIE}/src/Tools/Flux                               &&  $(MAKE) install && \
 	cd ${GENIE}/src/Tools/Geometry                           &&  $(MAKE) install && \
 	cd ${GENIE}/src/Tools/Masterclass                        &&  $(MAKE) install && \
+	cd ${GENIE}/src/Tools/VMC                                &&  $(MAKE) install && \
 	cd ${GENIE}
 
 
@@ -479,6 +494,7 @@ purge: FORCE
 	cd ${GENIE}/src/Tools/Flux                               &&  $(MAKE) purge && \
 	cd ${GENIE}/src/Tools/Geometry                           &&  $(MAKE) purge && \
 	cd ${GENIE}/src/Tools/Masterclass                        &&  $(MAKE) purge && \
+	cd ${GENIE}/src/Tools/VMC                                &&  $(MAKE) purge && \
 	cd ${GENIE}
 
 clean: clean-files clean-dir clean-etc
@@ -535,6 +551,7 @@ clean-files: FORCE
 	cd ${GENIE}/src/Tools/Flux                               &&  $(MAKE) clean && \
 	cd ${GENIE}/src/Tools/Geometry                           &&  $(MAKE) clean && \
 	cd ${GENIE}/src/Tools/Masterclass                        &&  $(MAKE) clean && \
+	cd ${GENIE}/src/Tools/VMC                                &&  $(MAKE) clean && \
 	cd ${GENIE}/src/Apps                                     &&  $(MAKE) clean && \
 	cd ${GENIE}/src/scripts                                  &&  $(MAKE) clean && \
 	cd ${GENIE}
@@ -604,6 +621,7 @@ distclean: FORCE
 	cd ${GENIE}/src/Tools/Flux                               &&  $(MAKE) distclean && \
 	cd ${GENIE}/src/Tools/Geometry                           &&  $(MAKE) distclean && \
 	cd ${GENIE}/src/Tools/Masterclass                        &&  $(MAKE) distclean && \
+	cd ${GENIE}/src/Tools/VMC                                &&  $(MAKE) distclean && \
 	cd ${GENIE}/src/Apps                                     &&  $(MAKE) distclean && \
 	cd ${GENIE}/src/scripts                                  &&  $(MAKE) distclean && \
 	cd ${GENIE}

config/EventGenerator.xml

@@ -115,6 +115,16 @@ XSecModel                   alg     Yes  Cross section model used at the thread
   </param_set>
 
 
+ <!--
+  VMC / EventLibraryInterface
+ -->
+  <param_set name="EventLibraryInterface">
+     <param type="string" name="VldContext"> </param>
+     <param type="int"    name="NModules">   1                                                    </param>
+     <param type="alg"    name="Module-0">   genie::vmc::EventLibraryInterface/Default            </param>
+     <param type="alg"    name="ILstGen">    genie::vmc::EvtLibInteractionListGenerator/Default   </param>
+  </param_set>
+
 
   <!--
   Event generation for Smith-Moniz CCQE model

config/EventGeneratorListAssembler.xml

@@ -47,6 +47,11 @@ Generator-%d                alg     No
 
 -->
 
+  <param_set name="EventLibraryInterface"> 
+     <param type="int" name="NGenerators">  1                                             </param>
+     <param type="alg" name="Generator-0">  genie::EventGenerator/EventLibraryInterface   </param>
+  </param_set>
+
   <param_set name="AnomalyMediatedNuGamma"> 
      <param type="int" name="NGenerators">  1                                  </param>
      <param type="alg" name="Generator-0">  genie::EventGenerator/AM-NUGAMMA   </param>

config/EventLibraryInterface.xml

@@ -0,0 +1,8 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<alg_conf>
+  <param_set name="Default">
+    <param type="string" name="LibraryPath"> Dummy.root </param>
+    <param type="bool" name="OnDemand"> true </param>
+  </param_set>
+</alg_conf>

config/EvtLibInteractionListGenerator.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<alg_conf>
+  <param_set name="Default">
+  </param_set>
+</alg_conf>

config/EvtLibPXSec.xml

@@ -0,0 +1,7 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+
+<alg_conf>
+  <param_set name="Default">
+    <param type="string" name="LibraryPath"> Dummy.root </param>
+  </param_set>
+</alg_conf>

config/G18_02a/ModelConfiguration.xml

@@ -138,7 +138,7 @@ STFC, Rutherford Appleton Laboratory
   <param type="alg" name="XSecModel@genie::EventGenerator/NucleonDecay">   genie::DummyPXSec/Default             </param>
   <param type="alg" name="XSecModel@genie::EventGenerator/NNBarOsc">       genie::NNBarOscDummyPXSec/Default     </param>
 
-
+  <param type="alg" name="XSecModel@genie::EventGenerator/EventLibraryInterface">  genie::vmc::EvtLibPXSec/Default   </param>
  </param_set>
 
 

config/master_config.xml

@@ -208,4 +208,8 @@
    <!-- ****** CONFIGURATION FOR REWEIGHT ALGORITHMS ****** -->
    <config alg="genie::rew::GSystUncertaintyTable">     GSystUncertaintyTable.xml        </config>
 
+ <!-- ****** VMC / EventLibraryInterface ****** -->
+   <config alg="genie::vmc::EventLibraryInterface">      EventLibraryInterface.xml       </config>
+   <config alg="genie::vmc::EvtLibInteractionListGenerator"> EvtLibInteractionListGenerator.xml  </config>
+   <config alg="genie::vmc::EvtLibPXSec"> EvtLibPXSec.xml  </config>
 </genie_config>

src/Tools/VMC/EventLibraryInterface.cxx

@@ -0,0 +1,261 @@
+//____________________________________________________________________________
+/*
+ Copyright (c) 2003-2020, The GENIE Collaboration
+ For the full text of the license visit http://copyright.genie-mc.org
+
+
+*/
+//____________________________________________________________________________
+
+#include "Framework/Messenger/Messenger.h"
+#include "Framework/Numerical/RandomGen.h"
+#include "Framework/GHEP/GHepRecord.h"
+#include "Framework/GHEP/GHepParticle.h"
+#include "Framework/ParticleData/PDGUtils.h"
+#include "Framework/ParticleData/PDGLibrary.h"
+#include "Framework/Interaction/Interaction.h"
+#include "Tools/VMC/EventLibraryInterface.h"
+#include "Tools/VMC/RecordList.h"
+#include "Tools/VMC/Utils.h"
+
+#include "TFile.h"
+
+using namespace genie;
+using namespace genie::vmc;
+
+//____________________________________________________________________________
+EventLibraryInterface::EventLibraryInterface() :
+  EventRecordVisitorI("genie::vmc::EventLibraryInterface"),
+  fRecordFile(0)
+{
+
+}
+//____________________________________________________________________________
+EventLibraryInterface::EventLibraryInterface(string config) :
+  EventRecordVisitorI("genie::vmc::EventLibraryInterface", config),
+  fRecordFile(0)
+{
+
+}
+//____________________________________________________________________________
+EventLibraryInterface::~EventLibraryInterface()
+{
+  for(auto it: fRecords) delete it.second;
+  delete fRecordFile;
+}
+//____________________________________________________________________________
+void EventLibraryInterface::ProcessEventRecord(GHepRecord * event) const
+{
+// Get event summary constructed by GENIE
+//
+  Interaction* interaction = event->Summary();
+  const InitialState & init_state = interaction->InitState();
+
+  const Record* rec = GetRecord(interaction);
+  if(!rec) return; // Reason has already been printed
+
+  std::unique_ptr<TLorentzVector> probe_p4(init_state.GetProbeP4(kRfLab));
+
+  unsigned int nLep = 0;
+  for(const Particle& part: rec->parts) if(pdg::IsLepton(part.pdg)) ++nLep;
+
+  const int firstLep = (nLep == 0) ? -1 : 2;
+  const int lastLep  = (nLep == 0) ? -1 : 1+nLep;
+
+  const int firstHad = (rec->parts.size() == nLep) ? -1 : 2+nLep;
+  const int lastHad  = (rec->parts.size() == nLep) ? -1 : 1+rec->parts.size();
+
+  // Neutrino is a parent to the lepton(s)
+  event->AddParticle(init_state.ProbePdg(),
+                     kIStInitialState,
+                     -1, -1, firstLep, lastLep,
+                     *probe_p4,
+                     TLorentzVector(0, 0, 0, 0));
+
+  const int tgt_A    = init_state.Tgt().A();
+  const int tgt_Z    = init_state.Tgt().Z();
+  const int tgt_pdgc = pdg::IonPdgCode(tgt_A, tgt_Z);
+  // REVIEW: it was like this in Costas' example code. Can it just be
+  // init_state.TgtPdg()?
+
+  LOG("ELI", pINFO)
+    << "Adding nucleus [A = " << tgt_A << ", Z = " << tgt_Z
+    << ", pdg = " << tgt_pdgc << "]";
+
+  // Nucleus is a parent to everything else
+  event->AddParticle(tgt_pdgc,
+                     kIStInitialState,
+                     -1, -1, firstHad, lastHad,
+                     0, 0, 0, PDGLibrary::Instance()->Find(tgt_pdgc)->Mass(),
+                     0, 0, 0, 0);
+
+  const std::vector<TVector3> basis = Basis(probe_p4->Vect());
+
+  for(bool lep: {true, false}){ // make sure lepton in first
+    for(const Particle& part: rec->parts){
+      if(lep != pdg::IsLepton(part.pdg)) continue;
+
+      // Fix up the outgoing lepton for NC events (in the library it's always
+      // nu_mu...)
+      int pdg = part.pdg;
+      if(interaction->ProcInfo().IsWeakNC()) pdg = init_state.ProbePdg();
+
+      event->AddParticle(pdg,
+                         kIStStableFinalState,
+                         (lep ? 0 : 1), -1, -1, -1, // child of the neutrino or nucleus
+                         TLorentzVector(part.px*basis[0] +
+                                        part.py*basis[1] +
+                                        part.pz*basis[2],
+                                        part.E),
+                         TLorentzVector(0, 0, 0, 0));
+    }
+  }
+}
+
+//____________________________________________________________________________
+const Record* EventLibraryInterface::GetRecord(const Interaction* interaction) const
+{
+  if(fRecords.empty()) LoadRecords();
+
+  const InitialState& init_state = interaction->InitState();
+
+  const double probe_E = init_state.ProbeE(kRfLab);
+
+  if(!init_state.Tgt().IsNucleus()){
+    LOG("ELI", pINFO) << "Skippping non-nuclear target " << init_state;
+    return 0;
+  }
+
+  const int tgt_pdgc = init_state.TgtPdg();
+
+  const ProcessInfo& proc = interaction->ProcInfo();
+
+  if(!proc.IsWeakCC() && !proc.IsWeakNC()){
+    LOG("ELI", pINFO) << "Skipping unknown process " << proc;
+    return 0;
+  }
+
+  int probe_pdgc = init_state.ProbePdg();
+
+  // Use nu_mu for NC by convention
+  if(proc.IsWeakNC()){
+    if(probe_pdgc > 0) probe_pdgc = +14; else probe_pdgc = -14;
+  }
+
+  const Key key(tgt_pdgc, probe_pdgc, proc.IsWeakCC());
+
+  const auto rec_it = fRecords.find(key);
+
+  if(rec_it == fRecords.end()){
+    LOG("ELI", pINFO) << "Skippping " << key << " -- not found in library";
+    return 0;
+  }
+
+  const Record* rec = rec_it->second->GetRecord(probe_E);
+
+  if(!rec){
+    LOG("ELI", pINFO) << "Skippping " << key << " at " << probe_E << " GeV -- not found in library";
+    return 0;
+  }
+
+  return rec;
+}
+
+//____________________________________________________________________________
+std::vector<TVector3> EventLibraryInterface::Basis(const TVector3& z) const
+{
+  TVector3 up(0, 1, 0);
+  if(up.Dot(z) == 0) up = TVector3(1, 0, 0);
+
+  const TVector3 x = up.Cross(z).Unit(); // Perpendicular to neutrino and up
+  const TVector3 y = x.Cross(z).Unit();  // Defines the third axis
+
+  const double a = RandomGen::Instance()->RndEvg().Uniform(0, 2*M_PI);
+
+  const TVector3 xp =  cos(a) * x + sin(a) * y;
+  const TVector3 yp = -sin(a) * x + cos(a) * y;
+  const TVector3 zp = z.Unit();
+
+  return {xp, yp, zp};
+}
+
+//____________________________________________________________________________
+void EventLibraryInterface::Configure(const Registry & config)
+{
+  Algorithm::Configure(config);
+}
+//___________________________________________________________________________
+void EventLibraryInterface::Configure(string config)
+{
+  Algorithm::Configure(config);
+}
+
+//___________________________________________________________________________
+void EventLibraryInterface::LoadRecords() const
+{
+  assert(!fRecordFile);
+
+  std::string libPath;
+  if(!GetParamDef( "LibraryPath", libPath, std::string() )){
+    LOG("ELI", pFATAL) << "Must specify 'LibraryPath'";
+    exit(1);
+  }
+
+  Expand(libPath);
+
+  bool onDemand;
+  GetParamDef( "OnDemand", onDemand, true );
+
+  PDGLibrary* pdglib = PDGLibrary::Instance();
+
+  fRecordFile = new TFile(libPath.c_str());
+  if(fRecordFile->IsZombie()) exit(1);
+
+  TIter next(fRecordFile->GetListOfKeys());
+  while(TObject* dir = next()){
+    const std::string& tgtName = dir->GetName();
+    const TParticlePDG* tgtPart = pdglib->DBase()->GetParticle(tgtName.c_str());
+    if(!tgtPart){
+      LOG("ELI", pWARN) << "Unknown nucleus " << tgtName
+                        << " found in " << libPath
+                        << " -- skipping";
+      continue;
+    }
+
+    for(int sign: {+1, -1}){
+      for(int pdg: {12, 14, 16}){
+        for(bool iscc: {true, false}){
+          // NCs should be the same for all flavours. Use numu by convention.
+          if(!iscc && pdg != 14) continue;
+
+          std::string nuName = pdglib->Find(sign*pdg)->GetName();
+          if(!iscc) nuName = (sign > 0) ? "nu" : "nu_bar";
+
+          const std::string treeName =
+            TString::Format("%s/%s/%s/records",
+                            tgtName.c_str(),
+                            iscc ? "cc" : "nc",
+                            nuName.c_str()).Data();
+
+          const Key key(tgtPart->PdgCode(), sign*pdg, iscc);
+
+          TTree* tr = (TTree*)fRecordFile->Get(treeName.c_str());
+
+          if(!tr){
+            LOG("ELI", pINFO) << treeName << " not found in "
+                              << libPath << " -- skipping";
+            continue;
+          }
+
+          if(onDemand)
+            fRecords[key] = new OnDemandRecordList(tr, treeName);
+          else
+            fRecords[key] = new SimpleRecordList(tr, treeName);
+        } // end for nucleus
+      } // end for pdg
+    } // end for sign
+  } // end for iscc
+
+  // Need to keep the record file open for OnDemand, but not Simple
+  if(!onDemand){delete fRecordFile; fRecordFile = 0;}
+}