Merge branch 'develop' into bugfix/ext

Author: kjplows

sbncode/CAFMaker/FillTrue.cxx

@@ -1260,38 +1260,47 @@ caf::Wall_t caf::GetWallCross(const geo::BoxBoundedGeo &volume, const TVector3 p
   TVector3 direction = (p1 - p0) * ( 1. / (p1 - p0).Mag());
   std::vector<TVector3> intersections = volume.GetIntersections(p0, direction);
 
-  assert(intersections.size() == 2);
+  /*
+   * There are either two, or one, or zero intersection points.
+   * As per larcorealg/Geometry/BoxBoundedGeo.h documentation:
+   * If the return std::vector is empty the trajectory does not intersect with the box.
+   * Normally the return value should have one (if the trajectory originates in the box) or two (else) entries.
+   * If the return value has two entries the first represents the entry point and the second the exit point
+   */
+
+  if( intersections.empty() ) return caf::kWallNone;
 
   // get the intersection point closer to p0
-  int intersection_i = ((intersections[0] - p0).Mag() < (intersections[1] - p0).Mag()) ? 0 : 1;
+  TVector3 closestIntersection;
+  double minDistance2 = std::numeric_limits<double>::max();
+
+  for( TVector3 const & point : intersections ) {
+    const double d2 = (point - p0).Mag2();
+    if( d2 > minDistance2 ) continue;
+    minDistance2 = d2;
+    closestIntersection = point;
+  }
 
   double eps = 1e-3;
-  if (abs(intersections[intersection_i].X() - volume.MinX()) < eps) {
-    //std::cout << "Left\n";
+  if (abs(closestIntersection.X() - volume.MinX()) < eps) {
     return caf::kWallLeft;
   }
-  else if (abs(intersections[intersection_i].X() - volume.MaxX()) < eps) {
-    //std::cout << "Right\n";
+  else if (abs(closestIntersection.X() - volume.MaxX()) < eps) {
     return caf::kWallRight;
   }
-  else if (abs(intersections[intersection_i].Y() - volume.MinY()) < eps) {
-    //std::cout << "Bottom\n";
+  else if (abs(closestIntersection.Y() - volume.MinY()) < eps) {
     return caf::kWallBottom;
   }
-  else if (abs(intersections[intersection_i].Y() - volume.MaxY()) < eps) {
-    //std::cout << "Top\n";
+  else if (abs(closestIntersection.Y() - volume.MaxY()) < eps) {
     return caf::kWallTop;
   }
-  else if (abs(intersections[intersection_i].Z() - volume.MinZ()) < eps) {
-    //std::cout << "Front\n";
+  else if (abs(closestIntersection.Z() - volume.MinZ()) < eps) {
     return caf::kWallFront;
   }
-  else if (abs(intersections[intersection_i].Z() - volume.MaxZ()) < eps) {
-    //std::cout << "Back\n";
+  else if (abs(closestIntersection.Z() - volume.MaxZ()) < eps) {
     return caf::kWallBack;
   }
   else assert(false);
-  //std::cout << "None\n";
 
   return caf::kWallNone;
 }//GetWallCross

sbncode/LArG4/CMakeLists.txt

@@ -1,15 +1,5 @@
 add_subdirectory(PhysicsLists)
 
-cet_build_plugin(G4InfoReducer art::module
-                       LIBRARIES
-                           larsim::Simulation
-                           larcore::Geometry_Geometry_service
-                           lardata::DetectorInfoServices_DetectorClocksServiceStandard_service
-                           art::Persistency_Common
-                           art::Persistency_Provenance
-                           art::Utilities canvas::canvas
-                           BASENAME_ONLY)
-
 set(
   MODULE_LIBRARIES
     art::Utilities
@@ -61,14 +51,5 @@ cet_build_plugin(MergeSimSourcesSBN art::EDProducer
   BASENAME_ONLY
 )
 
-cet_build_plugin(SimpleMerge art::EDProducer
-  LIBRARIES
-  lardataobj::MCBase
-  larcorealg::headers
-  larcorealg::CoreUtils
-  messagefacility::MF_MessageLogger
-)
-
 install_headers()
-install_fhicl()
 install_source()

sbncode/LArG4/g4inforeducer.fcl

@@ -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,7 +295,28 @@ 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;
+
+      // If Dk2Nu flux, use ancestors to evaluate tptype
+      if (fluxlist[inu].fFluxType == simb::kDk2Nu) {
+
+        for( const bsim::Ancestor & ancestor : dk2nu_v->at(inu).ancestor ) {
+          std::string aproc = ancestor.proc;
+          if( (aproc.find("BooNEHadronInelastic:BooNEpBeInteraction") != std::string::npos) && (aproc.find("QEBooNE") == std::string::npos) ) {
+            flux.ftptype = ancestor.pdg;
+            flux.ftpx = ancestor.startpx;
+            flux.ftpy = ancestor.startpy;
+            flux.ftpz = ancestor.startpz;
+          } // found it
+        } // find first inelastic
+      }
+
+      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);
           return weights;//done, all 1
@@ -301,19 +330,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/LArG4/simplemerge.fcl

@@ -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"

sbncode/SBNEventWeight/Calculators/BNBFlux/FluxCalcPrep.cxx

@@ -297,6 +297,9 @@ std::vector<float> Geant4WeightCalc::GetWeight(art::Event& e, size_t itruth ) {
         double Z = p.Position(i).Z();
         geo::Point_t testpoint1 { X, Y, Z };
         const TGeoMaterial* testmaterial1 = fGeometryService->Material( testpoint1 );
+	if( ! testmaterial1 ) {
+	  break; // The trajectory point is outside the world, or the geometry service is unable to handle the point
+	}
         //For now, just going to reweight the points within the LAr of the TPC
         // TODO check if this is right
         if ( !strcmp( testmaterial1->GetName(), "LAr" ) ){