Update GENIEGen and other ppfx-related updates

CMakeLists.txt

@@ -40,6 +40,10 @@ find_package(nurandom REQUIRED EXPORT)
 find_package(nusimdata REQUIRED EXPORT)
 find_package(nutools REQUIRED EXPORT)
 
+find_package(dk2nugenie REQUIRED EXPORT)
+find_package(dk2nudata REQUIRED EXPORT)
+find_package(ppfx REQUIRED EXPORT)
+
 find_package(Boost REQUIRED COMPONENTS headers math_tr1 REQUIRED EXPORT) # boost::math::policies, boost::math::ellint*
 find_package(CLHEP COMPONENTS Evaluator Geometry Random Vector REQUIRED EXPORT)
 find_package(CRY REQUIRED EXPORT) # find module in nutools

larsim/CMakeLists.txt

@@ -18,3 +18,4 @@ add_subdirectory(MergeSimSources)
 add_subdirectory(SimFilters)
 add_subdirectory(TriggerAlgo)
 add_subdirectory(PhotonHitConverter)
+add_subdirectory(PPFXFluxReader)

larsim/EventGenerator/GENIE/CMakeLists.txt

@@ -1,6 +1,7 @@
 cet_build_plugin(GENIEGen art::EDProducer
   LIBRARIES PRIVATE
   larcore::Geometry_Geometry_service
+  lardata::Utilities
   larcoreobj::SummaryData
   lardataalg::MCDumpers
   lardataobj::Simulation
@@ -14,6 +15,13 @@ cet_build_plugin(GENIEGen art::EDProducer
   canvas::canvas
   messagefacility::MF_MessageLogger
   fhiclcpp::fhiclcpp
+  GENIE::GFwEG
+  GENIE::GFwAlg
+  GENIE::GFwMsg
+  GENIE::GFwGHEP
+  GENIE::GTlFlx
+  dk2nu::Genie
+  dk2nu::Tree
   ROOT::Core
   ROOT::EG
   ROOT::Hist

larsim/EventGenerator/GENIE/GENIEGen_module.cc

@@ -38,13 +38,23 @@
 #include "larcore/Geometry/Geometry.h"
 #include "larcoreobj/SummaryData/POTSummary.h"
 #include "larcoreobj/SummaryData/RunData.h"
+#include "lardata/Utilities/AssociationUtil.h"
 #include "lardataalg/MCDumpers/MCDumpers.h" // sim::dump namespace
 #include "lardataobj/Simulation/BeamGateInfo.h"
 #include "nugen/EventGeneratorBase/GENIE/GENIEHelper.h"
 #include "nusimdata/SimulationBase/GTruth.h"
 #include "nusimdata/SimulationBase/MCFlux.h"
 #include "nusimdata/SimulationBase/MCTruth.h"
 
+// dk2nu extensions
+#include "dk2nu/genie/GDk2NuFlux.h"
+#include "dk2nu/tree/NuChoice.h"
+#include "dk2nu/tree/dk2nu.h"
+
+#include "GENIE/Framework/EventGen/EventRecord.h"
+#include "nugen/EventGeneratorBase/GENIE/EVGBAssociationUtil.h"
+#include "nugen/EventGeneratorBase/evgenbase.h"
+
 ///Event Generation using GENIE, cosmics or single particles
 namespace evgen {
   /**
@@ -166,6 +176,14 @@ namespace evgen {
     produces<art::Assns<simb::MCTruth, simb::GTruth>>();
     produces<std::vector<sim::BeamGateInfo>>();
 
+    // dk2nu additions
+    if (pset.get<std::string>("FluxType").find("dk2nu") != std::string::npos) {
+      produces<std::vector<bsim::Dk2Nu>>();
+      produces<std::vector<bsim::NuChoice>>();
+      produces<art::Assns<simb::MCTruth, bsim::Dk2Nu>>();
+      produces<art::Assns<simb::MCTruth, bsim::NuChoice>>();
+    }
+
     std::string beam_type_name = pset.get<std::string>("BeamName");
 
     if (beam_type_name == "numi")
@@ -374,6 +392,15 @@ namespace evgen {
     std::unique_ptr<std::vector<sim::BeamGateInfo>> gateCollection(
       new std::vector<sim::BeamGateInfo>);
 
+    std::unique_ptr<std::vector<bsim::Dk2Nu>> dk2nucol(new std::vector<bsim::Dk2Nu>);
+    std::unique_ptr<std::vector<bsim::NuChoice>> nuchoicecol(new std::vector<bsim::NuChoice>);
+    std::unique_ptr<art::Assns<simb::MCTruth, bsim::Dk2Nu>> dk2nuassn(
+      new art::Assns<simb::MCTruth, bsim::Dk2Nu>);
+    std::unique_ptr<art::Assns<simb::MCTruth, bsim::NuChoice>> nuchoiceassn(
+      new art::Assns<simb::MCTruth, bsim::NuChoice>);
+
+    genie::flux::GDk2NuFlux* dk2nuDriver =
+      dynamic_cast<genie::flux::GDk2NuFlux*>(fGENIEHelp->GetFluxDriver(true));
     while (truthcol->size() < 1) {
       while (!fGENIEHelp->Stop()) {
 
@@ -396,6 +423,21 @@ namespace evgen {
 
           FillHistograms(truth);
 
+          if (dk2nuDriver) {
+            const bsim::Dk2Nu& dk2nuObj = dk2nuDriver->GetDk2Nu();
+            dk2nucol->push_back(dk2nuObj);
+            const bsim::NuChoice& nuchoiceObj = dk2nuDriver->GetNuChoice();
+            nuchoicecol->push_back(nuchoiceObj);
+            util::CreateAssn(
+              evt, *truthcol, *dk2nucol, *dk2nuassn, dk2nucol->size() - 1, dk2nucol->size());
+            util::CreateAssn(evt,
+                             *truthcol,
+                             *nuchoicecol,
+                             *nuchoiceassn,
+                             nuchoicecol->size() - 1,
+                             nuchoicecol->size());
+          }
+
           // check that the process code is not unsupported by GENIE
           // (see issue #18025 for reference);
           // if it is, print all the information we can about this truth record
@@ -439,6 +481,14 @@ namespace evgen {
     evt.put(std::move(tgtassn));
     evt.put(std::move(gateCollection));
 
+    // dk2nu additions
+    if (dk2nuDriver) {
+      evt.put(std::move(dk2nucol));
+      evt.put(std::move(nuchoicecol));
+      evt.put(std::move(dk2nuassn));
+      evt.put(std::move(nuchoiceassn));
+    }
+
     return;
   }
 

larsim/EventWeight/Calculators/CMakeLists.txt

@@ -1,6 +1,21 @@
+
 cet_make_library(SOURCE GenieWeightCalc.cxx
+                        PPFXCVWeightCalc.cxx
+                        PPFXMIPPKaonWeightCalc.cxx
+                        PPFXMIPPPionWeightCalc.cxx
+                        PPFXOtherWeightCalc.cxx
+                        PPFXTargAttenWeightCalc.cxx
+                        PPFXThinKaonWeightCalc.cxx
+                        PPFXThinMesonWeightCalc.cxx
+                        PPFXThinNeutronPionWeightCalc.cxx
+                        PPFXThinNucAWeightCalc.cxx
+                        PPFXThinNucWeightCalc.cxx
+                        PPFXThinPionWeightCalc.cxx
+                        PPFXTotAbsorpWeightCalc.cxx
+                        PPFXWeightCalc.cxx
   LIBRARIES
   PRIVATE
+  ppfx::ppfx
   larsim::EventWeight_Base
   nugen::EventGeneratorBase_GENIE
   nusimdata::SimulationBase
@@ -9,6 +24,8 @@ cet_make_library(SOURCE GenieWeightCalc.cxx
   fhiclcpp::fhiclcpp
   cetlib_except::cetlib_except
   ${GENIE_LIB_LIST}
+  dk2nu::Tree
+  dk2nu::Genie
   CLHEP::Random
   log4cpp::log4cpp # FIXME Should be transitive from GENIE target(s).
 )

larsim/EventWeight/Calculators/PPFXCVWeightCalc.cxx

@@ -0,0 +1,166 @@
+
+#include "larsim/EventWeight/Base/WeightCalc.h"
+#include "larsim/EventWeight/Base/WeightCalcCreator.h"
+
+#include "CLHEP/Random/RandGaussQ.h"
+
+#include "MakeReweight.h"
+#include "cetlib/search_path.h"
+#include "dk2nu/tree/dk2nu.h"
+#include "dk2nu/tree/dkmeta.h"
+#include "nugen/EventGeneratorBase/GENIE/MCTruthAndFriendsItr.h"
+
+#include "TSystem.h"
+
+namespace evwgh {
+  class PPFXCVWeightCalc : public WeightCalc {
+  public:
+    PPFXCVWeightCalc();
+    void Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine) override;
+    std::vector<std::vector<double>> GetWeight(art::Event& e) override;
+
+  private:
+    std::string fGenieModuleLabel;
+
+    std::vector<std::string> fInputLabels;
+    std::string fPPFXMode;
+    std::string fMode;
+    std::string fHorn;
+    std::string fTarget;
+    int fSeed;
+    int fVerbose;
+    NeutrinoFluxReweight::MakeReweight* fPPFXrw;
+
+    DECLARE_WEIGHTCALC(PPFXCVWeightCalc)
+  };
+
+  PPFXCVWeightCalc::PPFXCVWeightCalc() {}
+
+  void PPFXCVWeightCalc::Configure(fhicl::ParameterSet const& p, CLHEP::HepRandomEngine& engine)
+  {
+    //get configuration for this function
+    fhicl::ParameterSet const& pset = p.get<fhicl::ParameterSet>(GetName());
+    fGenieModuleLabel = p.get<std::string>("genie_module_label");
+
+    //ppfx setup
+    fInputLabels = pset.get<std::vector<std::string>>("input_labels");
+    fPPFXMode = pset.get<std::string>("ppfx_mode");
+    fVerbose = pset.get<int>("verbose");
+    fMode = pset.get<std::string>("mode");
+    fHorn = pset.get<std::string>("horn_curr");
+    fTarget = pset.get<std::string>("target_config");
+    fSeed = pset.get<int>("random_seed", -1);
+
+    gSystem->Setenv("MODE", fPPFXMode.c_str());
+
+    fPPFXrw = NeutrinoFluxReweight::MakeReweight::getInstance();
+
+    std::string inputOptions;                               // Full path.
+    std::string mode_file = "inputs_" + fPPFXMode + ".xml"; // Just file name.
+    cet::search_path sp("FW_SEARCH_PATH");
+    sp.find_file(mode_file, inputOptions);
+
+    if (fSeed != -1) fPPFXrw->setBaseSeed(fSeed); // Set the random seed
+    std::cout << "is PPFX setup : " << fPPFXrw->AlreadyInitialized() << std::endl;
+    std::cout << "Setting PPFX, inputs: " << inputOptions << std::endl;
+    std::cout << "Setting Horn Current Configuration to: " << fHorn << std::endl;
+    std::cout << "Setting Target Configuration to: " << fTarget << std::endl;
+    if (!(fPPFXrw->AlreadyInitialized())) { fPPFXrw->SetOptions(inputOptions); }
+    std::cout << "PPFX just set with mode: " << fPPFXMode << std::endl;
+  }
+
+  std::vector<std::vector<double>> PPFXCVWeightCalc::GetWeight(art::Event& e)
+  {
+    std::vector<std::vector<double>> weight;
+    evgb::MCTruthAndFriendsItr mcitr(e, fInputLabels);
+    //calculate weight(s) here
+
+    int nmctruth = 0, nmatched = 0;
+    bool flag = true;
+    int ievt = -1;
+    std::vector<art::Handle<std::vector<bsim::Dk2Nu>>> dk2nus2 =
+      e.getMany<std::vector<bsim::Dk2Nu>>();
+    for (size_t dk2 = 0; dk2 < dk2nus2.size(); ++dk2) {
+      art::Handle<std::vector<bsim::Dk2Nu>> dk2nus = dk2nus2[dk2];
+    }
+    while ((flag = mcitr.Next())) {
+      std::string label = mcitr.GetLabel();
+      const simb::MCTruth* pmctruth = mcitr.GetMCTruth();
+      // const simb::GTruth*  pgtruth  = mcitr.GetGTruth();
+      //not-used//const simb::MCFlux*      pmcflux   = mcitr.GetMCFlux();
+      const bsim::Dk2Nu* pdk2nu = mcitr.GetDk2Nu();
+      //not-used//const bsim::NuChoice*    pnuchoice = mcitr.GetNuChoice();
+      // art::Ptr<simb::MCTruth>  mctruthptr = mcitr.GetMCTruthPtr();
+
+      ++ievt;
+      ++nmctruth;
+      if (fVerbose > 0) {
+        std::cout << "FluxWeightCalculator [" << std::setw(4) << ievt << "] "
+                  << " label \"" << label << "\" MCTruth@ " << pmctruth << " Dk2Nu@ " << pdk2nu
+                  << std::endl;
+      }
+      if (!pdk2nu) continue;
+      ++nmatched;
+
+      //RWH//bsim::Dk2Nu*  tmp_dk2nu  = fTmpDK2NUConversorAlg.construct_toy_dk2nu( &mctruth,&mcflux);
+      //RWH//bsim::DkMeta* tmp_dkmeta = fTmpDK2NUConversorAlg.construct_toy_dkmeta(&mctruth,&mcflux);
+      //RWH// those appear to have been memory leaks
+      //RWH// herein begins the replacment for the "construct_toy_dkmeta"
+
+      static bsim::DkMeta
+        dkmeta_obj; //RWH// create this on stack (destroyed when out-of-scope)  ... or static
+      dkmeta_obj.tgtcfg = fTarget;
+      dkmeta_obj.horncfg = fHorn;
+      (dkmeta_obj.vintnames).push_back("Index_Tar_In_Ancestry");
+      (dkmeta_obj.vintnames).push_back("Playlist");
+      bsim::DkMeta* tmp_dkmeta = &dkmeta_obj;
+      //RWH// herein ends this block
+
+      // sigh ....
+      //RWH// this is the signature in NeutrinoFluxReweight::MakeReweight :
+      //      //! create an interaction chain from the new dk2nu(dkmeta) format
+      //      void calculateWeights(bsim::Dk2Nu* nu, bsim::DkMeta* meta);
+      //RWH// these _should_ be "const <class>*" because we don't need to change them
+      //RWH// and the pointers we get out of the ART record are going to be const.
+      bsim::Dk2Nu* tmp_dk2nu = const_cast<bsim::Dk2Nu*>(pdk2nu); // remove const-ness
+      try {
+        fPPFXrw->calculateWeights(tmp_dk2nu, tmp_dkmeta);
+      }
+      catch (...) {
+        std::cerr << "Failed to calcualte weight" << std::endl;
+        continue;
+      }
+      //Get weights:
+      if (fMode == "reweight") {
+        double ppfx_cv_wgt = fPPFXrw->GetCVWeight();
+        std::vector<double> wvec = {ppfx_cv_wgt};
+        weight.push_back(wvec);
+      }
+      else {
+        std::vector<double> vmipppion = fPPFXrw->GetWeights("MIPPNumiPionYields");
+        std::vector<double> vmippkaon = fPPFXrw->GetWeights("MIPPNumiKaonYields");
+        std::vector<double> vtgtatt = fPPFXrw->GetWeights("TargetAttenuation");
+        std::vector<double> vabsorp = fPPFXrw->GetWeights("TotalAbsorption");
+        std::vector<double> vttpcpion = fPPFXrw->GetWeights("ThinTargetpCPion");
+        std::vector<double> vttpckaon = fPPFXrw->GetWeights("ThinTargetpCKaon");
+        std::vector<double> vttpcnucleon = fPPFXrw->GetWeights("ThinTargetpCNucleon");
+        std::vector<double> vttncpion = fPPFXrw->GetWeights("ThinTargetnCPion");
+        std::vector<double> vttnucleona = fPPFXrw->GetWeights("ThinTargetnucleonA");
+        std::vector<double> vttmesoninc = fPPFXrw->GetWeights("ThinTargetMesonIncident");
+        std::vector<double> vothers = fPPFXrw->GetWeights("Other");
+
+        std::vector<double> tmp_vhptot;
+        for (unsigned int iuniv = 0; iuniv < vtgtatt.size(); iuniv++) {
+          tmp_vhptot.push_back(float(vmipppion[iuniv] * vmippkaon[iuniv] * vtgtatt[iuniv] *
+                                     vabsorp[iuniv] * vttpcpion[iuniv] * vttpckaon[iuniv] *
+                                     vttpcnucleon[iuniv] * vttncpion[iuniv] * vttnucleona[iuniv] *
+                                     vttmesoninc[iuniv] * vothers[iuniv]));
+        }
+        weight.push_back(tmp_vhptot);
+      }
+    }
+
+    return weight;
+  }
+  REGISTER_WEIGHTCALC(PPFXCVWeightCalc)
+}