MiniBooNE Event Rescaling

sbncode/EventGenerator/MeVPrtl/Tools/Constants.cpp

@@ -22,7 +22,7 @@ Constants::Constants() {
   eta_mass = 0.547862; // GeV
   rho_mass = 0.77526; // GeV https://pdg.lbl.gov/2019/listings/rpp2019-list-rho-770.pdf
   etap_mass = 0.95778; // GeV
-
+  gamma_mass = 0;
   // Couplings
   fine_structure_constant = 7.2973525693e-3;// https://pdg.lbl.gov/2019/reviews/rpp2019-rev-phys-constants.pdf
   Gfermi = 1.166379e-5; // 1/GeV^2 https://pdg.lbl.gov/2020/reviews/rpp2020-rev-phys-constants.pdf
@@ -269,6 +269,8 @@ double PDG2Mass(int pdg) {
       return Constants::Instance().klong_mass;
     case 221:
       return Constants::Instance().eta_mass;
+    case 22:
+      return Constants::Instance().gamma_mass;
     case 11:
     case -11:
       return Constants::Instance().elec_mass;
@@ -280,6 +282,8 @@ double PDG2Mass(int pdg) {
     case 15:
     case -15:
       return Constants::Instance().tau_mass;
+    case 111:
+      return Constants::Instance().pizero_mass;
     case 211:
     case -211:
       return Constants::Instance().piplus_mass;

sbncode/EventGenerator/MeVPrtl/Tools/Constants.h

@@ -21,6 +21,7 @@ class Constants {
 
 public:
   // Masses
+  double gamma_mass;
   double elec_mass;
   double muon_mass;
   double piplus_mass;

sbncode/FluxReader/CMakeLists.txt

@@ -74,6 +74,7 @@ cet_make_library(
 cet_build_plugin(FluxReaderAna art::module LIBRARIES
                         sbncode_FluxReader
                         sbncode_FluxReader_FluxInterface
+			sbnobj::Common_SBNEventWeight
                         nusimdata::SimulationBase
                         larcoreobj::SummaryData
                         larcore::Geometry_Geometry_service

sbncode/FluxReader/FluxReaderAna_module.cc

@@ -35,6 +35,20 @@
 #include "nusimdata/SimulationBase/MCFlux.h"
 #include "nusimdata/SimulationBase/MCTruth.h"
 
+
+#include "canvas/Persistency/Common/FindManyP.h"
+#include "canvas/Persistency/Common/Ptr.h"
+#include "canvas/Persistency/Common/PtrVector.h"
+#include "sbnobj/Common/SBNEventWeight/EventWeightMap.h"
+
+
+#include <map>
+#include <vector>
+#include <string>
+#include <algorithm>
+#include <functional>
+#include <cmath>
+
 class FluxReaderAna;
 
 
@@ -101,6 +115,19 @@ class FluxReaderAna : public art::EDAnalyzer {
   float _nu_other_y; /// Y poisition of neutrino at the front face of the TPC (other)
   float _nu_other_z; /// Z poisition of neutrino at the front face of the TPC (other)
   float _nu_other_t; /// Time of the neutrino (other)
+
+
+  //For weighting:
+  std::string _flux_eventweight_multisim_producer;
+  bool _save_systs_flux;
+
+  int _evtwgt_flux_nfunc; /// Number of flux weight functions
+  std::vector<std::string> _evtwgt_flux_funcname; /// Names of weight functions
+  std::vector<int> _evtwgt_flux_nweight; /// Number of weights per function
+  std::vector<std::vector<float>> _evtwgt_flux_weight; /// Raw weights by function
+  std::vector<float> _evtwgt_flux_oneweight; /// Combined flux multisim weight
+
+
 };
 
 
@@ -122,6 +149,10 @@ FluxReaderAna::FluxReaderAna(fhicl::ParameterSet const& p)
   _x_cut = p.get<float>("XCut"); // cm
   _y_cut = p.get<float>("YCut"); // cm
 
+  _flux_eventweight_multisim_producer =
+    p.get<std::string>("FluxEventWeightProducer", "fluxweight");
+
+  _save_systs_flux = p.get<bool>("SaveSystsFlux", false);
 
   art::ServiceHandle<art::TFileService> fs;
   _tree = fs->make<TTree>("tree", "");
@@ -154,6 +185,14 @@ FluxReaderAna::FluxReaderAna(fhicl::ParameterSet const& p)
   _tree->Branch("nu_other_y", &_nu_other_y, "nu_other_y/F");
   _tree->Branch("nu_other_z", &_nu_other_z, "nu_other_z/F");
   _tree->Branch("nu_other_t", &_nu_other_t, "nu_other_t/F");
+
+
+  _tree->Branch("evtwgt_flux_nfunc", &_evtwgt_flux_nfunc, "evtwgt_flux_nfunc/I");
+  _tree->Branch("evtwgt_flux_funcname", "std::vector<std::string>", &_evtwgt_flux_funcname);
+  _tree->Branch("evtwgt_flux_nweight", "std::vector<int>", &_evtwgt_flux_nweight);
+  _tree->Branch("evtwgt_flux_weight", "std::vector<std::vector<float>>", &_evtwgt_flux_weight);
+  _tree->Branch("evtwgt_flux_oneweight", "std::vector<float>", &_evtwgt_flux_oneweight);
+
 }
 
 void FluxReaderAna::analyze(art::Event const& e)
@@ -167,6 +206,12 @@ void FluxReaderAna::analyze(art::Event const& e)
   e.getByLabel(_flux_label, mctruthHandle);
   std::vector<simb::MCTruth> const& mclist = *mctruthHandle;
 
+  std::vector<art::Ptr<simb::MCTruth>> mct_v;
+  art::fill_ptr_vector(mct_v, mctruthHandle);
+
+  art::FindManyP<sbn::evwgh::EventWeightMap> mct_to_fluxewm(mctruthHandle, e, 
+							    _flux_eventweight_multisim_producer);
+
   for(unsigned int inu = 0; inu < mclist.size(); inu++) {
     simb::MCParticle nu = mclist[inu].GetNeutrino().Nu();
     simb::MCFlux flux = fluxlist[inu];
@@ -228,6 +273,76 @@ void FluxReaderAna::analyze(art::Event const& e)
       }
     }
 
+    // Reset flux systematic containers every neutrino
+    _evtwgt_flux_nfunc = 0;
+    _evtwgt_flux_funcname.clear();
+    _evtwgt_flux_nweight.clear();
+    _evtwgt_flux_weight.clear();
+    _evtwgt_flux_oneweight.clear();
+
+    if (_save_systs_flux) {
+
+      std::vector<art::Ptr<sbn::evwgh::EventWeightMap>> flux_ewm_v = mct_to_fluxewm.at(inu);
+
+      if (flux_ewm_v.size() == 1) {
+
+	std::map<std::string, std::vector<float>> evtwgt_map = *(flux_ewm_v[0]);
+
+	std::vector<float> previous_weights;
+	std::vector<float> final_weights;
+
+        int countFunc = 0;
+
+        for (auto const& it : evtwgt_map) {
+	  std::string const& func_name = it.first;
+	  std::vector<float> const& weight_v = it.second;
+
+          if (previous_weights.empty()) {
+            previous_weights.resize(weight_v.size(), 1.f);
+            final_weights.resize(weight_v.size(), 1.f);
+          }
+
+          // Protect against mismatched universe counts
+          if (weight_v.size() != previous_weights.size()) {
+	    mf::LogWarning("FluxReaderAna")
+              << "Weight function " << func_name
+              << " has size " << weight_v.size()
+              << " but previous combined weight size is "
+              << previous_weights.size()
+              << ". Skipping this function.";
+            continue;
+          }
+
+          _evtwgt_flux_funcname.push_back(func_name);
+          _evtwgt_flux_weight.push_back(weight_v);
+          _evtwgt_flux_nweight.push_back((int)weight_v.size());
+          countFunc++;
+
+	  std::transform(previous_weights.begin(), previous_weights.end(),
+                         weight_v.begin(),
+                         final_weights.begin(),
+                         std::multiplies<float>());
+
+          previous_weights = final_weights;
+        }
+
+        _evtwgt_flux_nfunc = countFunc;
+        _evtwgt_flux_oneweight = final_weights;
+      }
+      else if (flux_ewm_v.empty()) {
+	mf::LogWarning("FluxReaderAna")
+          << "No associated EventWeightMap found for neutrino index "
+          << inu << " from producer " << _flux_eventweight_multisim_producer;
+      }
+      else {
+	mf::LogWarning("FluxReaderAna")
+          << "Found " << flux_ewm_v.size()
+          << " associated EventWeightMaps for neutrino index " << inu
+          << " from producer " << _flux_eventweight_multisim_producer
+          << ", expected 1.";
+      }
+    }
+
     _tree->Fill();
   }
 }

sbncode/FluxReader/job/run_fluxreader_dk2nu_sbnd.fcl

@@ -0,0 +1,7 @@
+# Driver fcl file for reading in BooNE files
+# for SBND
+#include "genie_sbnd.fcl"
+#include "run_fluxreader_sbnd.fcl"
+
+source.inputType: "dk2nu"
+

sbncode/FluxReader/job/run_fluxreader_sbnd.fcl

@@ -23,11 +23,68 @@ source:
   module_type:	    FluxReader
   skipEvents:       0
   maxEvents:        -1
-  inputType:        "gsimple"
+  inputType:        "dk2nu"
+  dk2nuConfig: "dk2nu_sbnd_v3"
   nBins:            200
   Elow:             0
   Ehigh:            10
   SelfIncrementRun: true
+  dk2nu_sbnd: {
+        rotmatrix: [ 1, 0, 0, 
+		     0, 1, 0, 
+		     0, 0, 1]
+        userbeam: [ 73.78, 0, 11000]
+
+        #x_detAV: [-200, 200]
+        #y_detAV: [-200, 200]
+        #z_detAV: [0, 500]
+
+	windowBase: [  500,  500, -1000 ]
+        window1:    [  500, -500, -1000 ]
+        window2:    [ -500,  500, -1000 ]
+
+	}	
+  dk2nu_sbnd_v2: {
+        rotmatrix: [ 1, 0, 0, 
+		     0, 1, 0, 
+		     0, 0, 1]
+        userbeam: [ 0, 0, 0, 73.78, 0, 11000]
+
+
+	windowBase: [  500,  500, -1000 ]
+        window1:    [  0, -1000, 0 ]
+        window2:    [ -1000,  0, 0 ]
+
+	}
+dk2nu_sbnd_v3: {
+         userbeam: [ 73.78, 0, 11000 ]
+         rotmatrix:  [ 1,  0,  0,
+                       0,  1,  0,
+                       0,  0,  1  ]
+         windowBase: [   500,  500, -1000 ]
+         window1:    [  -500,  500, -1000 ]
+         window2:    [   500, -500, -1000 ]
+         }
+
+dk2nu_sbnd_leo: {
+         userbeam: [ -73.78, 0, 11000 ]
+         rotmatrix:  [ 1,  0,  0,
+                       0,  1,  0,
+                       0,  0,  1  ]
+         windowBase: [   500,  500, -1000 ]
+         window1:    [  -500,  500, -1000 ]
+         window2:    [   500, -500, -1000 ]
+         }
+  dk2nu_bnb_at_uboone_v1: {
+        userbeam:   [ -130., 0., 47000 ]
+	#userbeam: [ -73.78, 0., -11000.]
+	rotmatrix:  [ 1,  0,  0,
+                   0,  1,  0,
+		   0,  0,  1  ]
+	windowBase: [  630,  500, -1000 ]
+	window1:    [  630, -500, -1000 ]
+	window2:    [ -370,  500, -1000 ]
+  }
 }
 
 outputs:
@@ -38,7 +95,7 @@ outputs:
     fileName:    "fluxreader.root"
     compressionLevel: 1
     dataTier:    "simulated"
-    SelectEvents: [ filter ]
+    SelectEvents: [ filter  ]
   }
 }
 

sbncode/FluxReader/job/run_fluxreader_sbnd_SBNEventWeight.fcl

@@ -38,6 +38,7 @@ physics.producers.fluxweight.weight_functions_flux: [horncurrent
                                                    , kzero
                                                    , piplus
                                                    , piminus
+						   , fluxhist
                                                      ]
 
 # Need to overwrite the parameter weight_functions as weight_functions_flux;

sbncode/FluxReader/job/run_fluxreader_sbnd_boone_SBNEventWeight.fcl

@@ -0,0 +1,50 @@
+# Driver fcl file for reading in gsimple files
+# for SBND
+# Add SBNEventWeight -- Keng Lin June 2021
+
+#include "seedservice.fcl"
+#include "services_sbnd.fcl"
+#include "eventweight_flux_sbn.fcl"
+#include "eventweight_genie_sbn.fcl"
+
+#include "run_fluxreader_sbnd.fcl"
+
+
+source.inputType: "boone"
+
+physics.producers.fluxweight:  @local::sbn_eventweight_flux
+# physics.producers.genieweight: @local::sbn_eventweight_genie
+
+physics.filter: [ rns
+                , fluxfilter
+                , fluxweight
+                # , genieweight
+                ]
+
+
+# physics.filters.fluxfilter.volumes: ["volWorld"]
+
+physics.producers.fluxweight.generator_module_label: flux
+
+# Customize what kind of weights we need, i.e. define weight_functions_flux:
+physics.producers.fluxweight.weight_functions_flux: [horncurrent
+                                                   , expskin
+                                                   , pioninexsec
+                                                   , pionqexsec
+                                                   , piontotxsec
+                                                   , nucleoninexsec
+                                                   , nucleonqexsec
+                                                   , nucleontotxsec
+                                                   , kplus
+                                                   , kminus
+                                                   , kzero
+                                                   , piplus
+                                                   , piminus
+						   , fluxhist
+                                                     ]
+
+# Need to overwrite the parameter weight_functions as weight_functions_flux;
+physics.producers.fluxweight.weight_functions: @local::physics.producers.fluxweight.weight_functions_flux
+
+
+

sbncode/FluxReader/job/run_fluxreader_sbnd_dk2nu_SBNEventWeight.fcl

@@ -0,0 +1,49 @@
+# Driver fcl file for reading in gsimple files
+# for SBND
+# Add SBNEventWeight -- Keng Lin June 2021
+
+#include "seedservice.fcl"
+#include "services_sbnd.fcl"
+#include "eventweight_flux_sbn.fcl"
+#include "eventweight_genie_sbn.fcl"
+
+#include "run_fluxreader_sbnd.fcl"
+
+source.inputType: "dk2nu"
+
+physics.producers.fluxweight:  @local::sbn_eventweight_flux
+# physics.producers.genieweight: @local::sbn_eventweight_genie
+
+physics.filter: [ rns
+                , fluxfilter
+                , fluxweight
+                # , genieweight
+                ]
+
+
+# physics.filters.fluxfilter.volumes: ["volWorld"]
+
+physics.producers.fluxweight.generator_module_label: flux
+
+# Customize what kind of weights we need, i.e. define weight_functions_flux:
+physics.producers.fluxweight.weight_functions_flux: [horncurrent
+                                                   , expskin
+                                                   , pioninexsec
+                                                   , pionqexsec
+                                                   , piontotxsec
+                                                   , nucleoninexsec
+                                                   , nucleonqexsec
+                                                   , nucleontotxsec
+                                                   , kplus
+                                                   , kminus
+                                                   , kzero
+                                                   , piplus
+                                                   , piminus
+						   , fluxhist
+                                                     ]
+
+# Need to overwrite the parameter weight_functions as weight_functions_flux;
+physics.producers.fluxweight.weight_functions: @local::physics.producers.fluxweight.weight_functions_flux
+
+
+