Make flux calculator compatible with g4bnb flux

Author: marcodeltutto

sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx

@@ -213,8 +213,16 @@ namespace sbn {
 //      std::cout<<"SBNEventWeight : getweight for the "<<inu<<" th particles of an event"<< std::endl;
       //MCFlux & MCTruth
       art::Handle< std::vector<simb::MCFlux> > mcFluxHandle;
-      e.getByLabel(fGeneratorModuleLabel,mcFluxHandle);
+      e.getByLabel(fGeneratorModuleLabel, mcFluxHandle);
       std::vector<simb::MCFlux> const& fluxlist = *mcFluxHandle;
+
+      art::Handle< std::vector<bsim::Dk2Nu> > dk2nuHandle;
+      std::vector<bsim::Dk2Nu> const* dk2nu_v = nullptr;
+      e.getByLabel(fGeneratorModuleLabel, dk2nuHandle);
+      if (dk2nuHandle.isValid()) {
+        dk2nu_v = dk2nuHandle.product(); 
+      }
+
       //or do the above 3 lines in one line
 //      auto const& mclist = *e.getValidHandle<std::vector<simb::MCTruth>>(fGeneratorModuleLabel);
 
@@ -287,9 +295,35 @@ namespace sbn {
 
         // First let's check that the parent of the neutrino we are looking for is 
         //  the particle we intended it to be, if not set all weights to 1
-      if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(fluxlist[inu].ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1
+
+      simb::MCFlux flux;
+      flux.ftptype = fluxlist[inu].ftptype;
+      flux.ftpx = fluxlist[inu].ftpx;
+      flux.ftpy = fluxlist[inu].ftpy;
+      flux.ftpz = fluxlist[inu].ftpz;
+
+      // int tptype = fluxlist[inu].ftptype;
+
+      // If Dk2Nu flux, use ancestors to evaluate tptype
+      if (fluxlist[inu].fFluxType == simb::kDk2Nu) {
+
+        // Find the first inelastic interaction
+        int firstInelastic = 0;
+        while (dk2nu_v->at(inu).ancestor[firstInelastic].proc.find("HadronInelastic")==std::string::npos||dk2nu_v->at(inu).ancestor[firstInelastic].proc.find("QEBooNE")!=std::string::npos) firstInelastic++;
+        
+        // Get the parent/ancestor (this should be the secondary in the p+Be interaction)        
+        flux.ftptype = dk2nu_v->at(inu).ancestor[firstInelastic].pdg;
+        flux.ftpx = dk2nu_v->at(inu).ancestor[firstInelastic].startpx;
+        flux.ftpy = dk2nu_v->at(inu).ancestor[firstInelastic].startpy;
+        flux.ftpz = dk2nu_v->at(inu).ancestor[firstInelastic].startpz;
+
+      }
+
+
+      if (std::find(fprimaryHad.begin(), fprimaryHad.end(),(flux.ftptype)) == fprimaryHad.end() ){//if it does not contain any particles we need get 1
           weights.resize( NUni);
           std::fill(weights.begin(), weights.end(), 1);
+          std::cout << "We don't need this parent, returning." << std::endl;
           return weights;//done, all 1
         }// Hadronic parent check
 
@@ -301,19 +335,18 @@ namespace sbn {
             std::vector< float > Vrandom = (fParameterSet.fParameterMap.begin())->second;//vector of random #
             std::vector< float > subVrandom;//sub-vector of random numbers;
             if( CalcType == "PrimaryHadronNormalization"){//Normalization
-              test_weight = PHNWeightCalc(fluxlist[inu], Vrandom[i]);
+              test_weight = PHNWeightCalc(flux, Vrandom[i]);
 
             } else {
               subVrandom = {Vrandom.begin()+i*FitCov->GetNcols(), Vrandom.begin()+(i+1)*FitCov->GetNcols()};
               if( CalcType == "PrimaryHadronFeynmanScaling"){//FeynmanScaling
-                test_weight = PHFSWeightCalc(fluxlist[inu], subVrandom);
+                test_weight = PHFSWeightCalc(flux, subVrandom);
 
               } else if( CalcType == "PrimaryHadronSanfordWang"){//SanfordWang
-                test_weight = PHSWWeightCalc(fluxlist[inu], subVrandom);
+                test_weight = PHSWWeightCalc(flux, subVrandom);
 
               } else if( CalcType == "PrimaryHadronSWCentralSplineVariation"){//SWCentaralSplineVariation
-                test_weight = PHSWCSVWeightCalc(fluxlist[inu], subVrandom);
-
+                test_weight = PHSWCSVWeightCalc(flux, subVrandom);
               } else throw cet::exception(__PRETTY_FUNCTION__) << GetName() << ": this shouldnt happen.."<<std::endl;
             }
 

sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.h

@@ -5,6 +5,10 @@
 #include "nusimdata/SimulationBase/MCFlux.h" 
 #include "nusimdata/SimulationBase/MCTruth.h"
 
+#include "dk2nu/tree/NuChoice.h"
+#include "dk2nu/tree/dk2nu.h"
+#include "dk2nu/tree/dkmeta.h"
+
 #include "sbncode/SBNEventWeight/Base/WeightCalc.h"
 #include "art/Framework/Principal/Event.h"
 #include "sbncode/SBNEventWeight/Base/WeightCalcCreator.h"