Feature/howard for numi2023a from gputnam configs

sbncode/CAFMaker/FillTrigger.cxx

@@ -13,6 +13,15 @@ namespace caf
     triggerInfo.global_trigger_det_time = trig.TriggerTime();
     double diff_ts = triggerInfo.global_trigger_det_time - triggerInfo.beam_gate_det_time;
     triggerInfo.trigger_within_gate = diff_ts;
+
+    triggerInfo.prev_global_trigger_time = addltrig_info.previousTriggerTimestamp;
+    triggerInfo.source_type = (int)addltrig_info.sourceType;
+    triggerInfo.trigger_type = (int)addltrig_info.triggerType;
+    triggerInfo.trigger_id = addltrig_info.triggerID;
+    triggerInfo.gate_id = addltrig_info.gateID;
+    triggerInfo.trigger_count = addltrig_info.triggerCount;
+    triggerInfo.gate_count = addltrig_info.gateCount;
+    triggerInfo.gate_delta = (int)addltrig_info.gateCountFromPreviousTrigger;
   }
 
   void FillTriggerMC(double absolute_time, caf::SRTrigger& triggerInfo) {

sbncode/CAFMaker/FillTrue.cxx

@@ -470,7 +470,9 @@ namespace caf {
     for(const caf::SRTrueParticle& part: srparticles){
       // save the G4 particles that came from this interaction
       if(part.interaction_id == (int)i) {
-        if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part);
+        //if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part);
+        // Feb. 15th 2024, Jaesung Kim: Saving non-primaries for G4Reweight study
+        srneutrino.prim.push_back(part);
 
         // total up the deposited energy
         for(int p = 0; p < 3; ++p) { 

sbncode/Calibration/TrackCaloSkimmer_module.cc

@@ -831,6 +831,41 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,
     }
   }
 
+  // Save the true trajectory
+  for (unsigned i_traj = 0; i_traj < particle.NumberTrajectoryPoints(); i_traj++) {
+    // Get trajectory point
+    TVector3 traj = particle.Position(i_traj).Vect();
+    geo::Point_t traj_p(traj.X(), traj.Y(), traj.Z());
+
+    // lookup TPC
+    geo::TPCID tpc; // invalid by default
+    for (auto const &cryo: geo->Iterate<geo::CryostatGeo>()) {
+      for (auto const& TPC : geo->Iterate<geo::TPCGeo>(cryo.ID())) {
+        if (TPC.ActiveBoundingBox().ContainsPosition(traj_p)) {
+          tpc = TPC.ID();
+          break;
+        }
+      }
+      if (tpc.isValid) break;
+    }
+
+    // add in space-charge-deflected position if applicable
+    geo::Point_t traj_p_sce = tpc.isValid ? TrajectoryToWirePosition(traj_p, tpc) : traj_p;
+
+    sbn::Vector3D traj_v;
+    traj_v.x = traj_p.x();
+    traj_v.y = traj_p.y();
+    traj_v.z = traj_p.z();
+
+    sbn::Vector3D traj_v_sce;
+    traj_v_sce.x = traj_p_sce.x();
+    traj_v_sce.y = traj_p_sce.y();
+    traj_v_sce.z = traj_p_sce.z();
+
+    trueparticle.traj.push_back(traj_v);
+    trueparticle.traj_sce.push_back(traj_v_sce);
+  }
+
   return trueparticle;
 }
 

sbncode/EventGenerator/MeVPrtl/Tools/ALP/Meson2ALP_tool.cc

@@ -76,6 +76,7 @@ class Meson2ALP : public IMeVPrtlFlux
   double fcG; //!< Axion coupling to gluon
   double fcB; //!< Axion coupling to U(1) B boson (before EW symmetry breaking)
   double fcW; //!< Axion coupling to SU(2) W bosons (before EW symmetry breaking)
+  bool fIncludeMassSuppressionFactor; // Whether to include a phenomonological suppression to mixing at high mass
 
   // branching ratios
   double fPi0BR;
@@ -115,6 +116,7 @@ void Meson2ALP::configure(fhicl::ParameterSet const &pset)
   fMaxWeightPi0 = pset.get<double>("MaxWeightPi0", 1.);
   fMaxWeightEta = pset.get<double>("MaxWeightEta", 1.);
   fMaxWeightEtaP = pset.get<double>("MaxWeightEtaP", 1.);
+  fIncludeMassSuppressionFactor = pset.get<bool>("IncludeMassSuppressionFactor", false);
 
   fPi0BR = Pi0BR();
   fEtaBR = EtaBR();
@@ -136,7 +138,7 @@ double Meson2ALP::EtaBR() const {
 
   double mixing_angle = (1. / sqrt(6)) * (fcG * fpion / ffa) * (fM*fM - (4./9.)*pzero_mass*pzero_mass) / (fM*fM - eta_mass*eta_mass);
   double qcd_rate_f = 1.;
-  if (fM > eta_mass) qcd_rate_f = pow(fM/eta_mass, -1.6);
+  if (fIncludeMassSuppressionFactor && fM > eta_mass) qcd_rate_f = pow(fM/eta_mass, -1.6);
 
   return mixing_angle*mixing_angle*qcd_rate_f;
 }
@@ -148,7 +150,7 @@ double Meson2ALP::EtaPBR() const {
 
   double mixing_angle = (1. / sqrt(12)) * (fcG * fpion / ffa) * (fM*fM - (16./9.)*pzero_mass*pzero_mass) / (fM*fM - etap_mass*etap_mass);
   double qcd_rate_f = 1.;
-  if (fM > etap_mass) qcd_rate_f = pow(fM/etap_mass, -1.6);
+  if (fIncludeMassSuppressionFactor && fM > etap_mass) qcd_rate_f = pow(fM/etap_mass, -1.6);
 
   return mixing_angle*mixing_angle*qcd_rate_f;
 }
@@ -160,7 +162,7 @@ double Meson2ALP::Pi0BR() const {
   double mixing_angle = (1./6) * (fcG * fpion / ffa) * fM*fM / (fM*fM - pzero_mass*pzero_mass);
 
   double qcd_rate_f = 1.;
-  if (fM > pzero_mass) qcd_rate_f = pow(fM/pzero_mass, -1.6);
+  if (fIncludeMassSuppressionFactor && fM > pzero_mass) qcd_rate_f = pow(fM/pzero_mass, -1.6);
 
   return mixing_angle*mixing_angle*qcd_rate_f;
 }
@@ -177,6 +179,10 @@ bool Meson2ALP::MakeFlux(const simb::MCFlux &flux, evgen::ldm::MeVPrtlFlux &alp,
 
   // Energy is same as for meson (don't worry about momentum conservation)
   double alp_energy = meson.mom.E();
+
+  // ignore if alp isn't energetic enough
+  if (alp_energy < fM) return false;
+
   double alp_momentum = sqrt(alp_energy*alp_energy - fM*fM);
 
   // Momentum 4-vector

sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.cxx

@@ -1,6 +1,9 @@
 #include "TrajectoryMCSFitter.h"
 #include "lardataobj/RecoBase/Track.h"
 #include "larcorealg/Geometry/geo_vectors_utils.h"
+#include "larcore/Geometry/Geometry.h"
+#include "larcore/CoreUtils/ServiceUtil.h"
+#include "larcorealg/Geometry/GeometryCore.h"
 #include "TMatrixDSym.h"
 #include "TMatrixDSymEigen.h"
 
@@ -15,7 +18,8 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
   vector<size_t> breakpoints;
   vector<float> segradlengths;
   vector<float> cumseglens;
-  breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens);
+  vector<bool> breakpointsgood;
+  breakTrajInSegments(traj, breakpoints, segradlengths, cumseglens, breakpointsgood);
   //
   // Fit segment directions, and get 3D angles between them
   //
@@ -28,7 +32,11 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
     if (p>0) {
       if (segradlengths[p]<-100. || segradlengths[p-1]<-100.) {
         dtheta.push_back(-999.);
-            } else { 
+      } 
+      if (!breakpointsgood[p] || !breakpointsgood[p-1]) {
+        dtheta.push_back(-999.);
+      }
+      else { 
         const double cosval = pcdir0.X()*pcdir1.X()+pcdir0.Y()*pcdir1.Y()+pcdir0.Z()*pcdir1.Z();
         //assert(std::abs(cosval)<=1);
         //units are mrad
@@ -49,15 +57,21 @@ recob::MCSFitResult TrajectoryMCSFitter::fitMcs(const recob::TrackTrajectory& tr
   }
   const ScanResult fwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenFwd, true,  momDepConst, pid);
   const ScanResult bwdResult = doLikelihoodScan(dtheta, segradlengths, cumLenBwd, false, momDepConst, pid);
-  //
+  // std::cout << "fwdResult.p=" << fwdResult.p << " fwdResult.pUnc=" << fwdResult.pUnc << " fwdResult.logL=" << fwdResult.logL << std::endl;
   return recob::MCSFitResult(pid,
                             fwdResult.p,fwdResult.pUnc,fwdResult.logL,
                             bwdResult.p,bwdResult.pUnc,bwdResult.logL,
                             segradlengths,dtheta);
 }
 
-void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector<size_t>& breakpoints, vector<float>& segradlengths, vector<float>& cumseglens) const {
-  //
+void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj, vector<size_t>& breakpoints, vector<float>& segradlengths, vector<float>& cumseglens, vector<bool>& breakpointsgood) const {
+  // set fiducial volume
+  std::vector<geo::BoxBoundedGeo> fiducialVolumes;
+  fiducialVolumes = setFiducialVolumes();
+  // set excluded volumes
+  std::vector<geo::BoxBoundedGeo> excludeVolumes;
+  excludeVolumes = setExcludeVolumes();
+
   const double trajlen = traj.Length();
   const int nseg = std::max(minNSegs_,int(trajlen/segLen_));
   const double thisSegLen = trajlen/double(nseg);
@@ -69,6 +83,8 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
   auto nextValid=traj.FirstValidPoint();
   breakpoints.push_back(nextValid);
   auto pos0 = traj.LocationAtPoint(nextValid);
+  bool pos0good = isInVolume(fiducialVolumes, pos0) && !isInVolume(excludeVolumes, pos0);
+  breakpointsgood.push_back(pos0good);
   nextValid = traj.NextValidPoint(nextValid+1);
   int npoints = 0;
   while (nextValid!=recob::TrackTrajectory::InvalidIndex) {
@@ -77,7 +93,9 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
     pos0=pos1;
     npoints++;
     if (thislen>=thisSegLen) {
+      bool pos1good = isInVolume(fiducialVolumes, pos1) && !isInVolume(excludeVolumes, pos1);
       breakpoints.push_back(nextValid);
+      breakpointsgood.push_back(pos1good);
       if (npoints>=minHitsPerSegment_) segradlengths.push_back(thislen*lar_radl_inv);
       else segradlengths.push_back(-999.);
       cumseglens.push_back(cumseglens.back()+thislen);
@@ -88,7 +106,10 @@ void TrajectoryMCSFitter::breakTrajInSegments(const recob::TrackTrajectory& traj
   }
   //then add last segment
   if (thislen>0.) {
+    auto endpointpos = traj.LocationAtPoint(traj.LastValidPoint());
+    bool endpointposgood = isInVolume(fiducialVolumes, endpointpos) && !isInVolume(excludeVolumes, endpointpos);
     breakpoints.push_back(traj.LastValidPoint()+1);
+    breakpointsgood.push_back(endpointposgood);
     segradlengths.push_back(thislen*lar_radl_inv);
     cumseglens.push_back(cumseglens.back()+thislen);
   }
@@ -202,7 +223,6 @@ double TrajectoryMCSFitter::mcsLikelihood(double p, double theta0x, std::vector<
   double result = 0;
   for (int i = beg; i != end; i+=incr ) {
     if (dthetaij[i]<0) {
-      //cout << "skip segment with too few points" << endl;
       continue;
     }
     //
@@ -302,3 +322,75 @@ double TrajectoryMCSFitter::GetE(const double initial_E, const double length_tra
   }
   return current_E;
 }
+
+std::vector<geo::BoxBoundedGeo> TrajectoryMCSFitter::setFiducialVolumes() const {
+  std::vector<geo::BoxBoundedGeo> fiducialVolumes;
+  std::vector<std::vector<geo::BoxBoundedGeo>> TPCVolumes;
+  std::vector<geo::BoxBoundedGeo> ActiveVolumes;
+
+  const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
+  for (auto const &cryoID: geometry->Iterate<geo::CryostatID>()) {
+    std::vector<geo::BoxBoundedGeo> thisTPCVolumes;
+    for (auto const& tpc: geometry->Iterate<geo::TPCGeo>(cryoID)) {
+      thisTPCVolumes.push_back(tpc.ActiveBoundingBox());
+    }
+    TPCVolumes.push_back(std::move(thisTPCVolumes));
+  }
+  for (const std::vector<geo::BoxBoundedGeo> &tpcs: TPCVolumes) {
+    double xMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinX() < rhs.MinX(); })->MinX();
+    double xMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxX() < rhs.MaxX(); })->MaxX();
+    double yMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinY() < rhs.MinY(); })->MinY();
+    double yMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxY() < rhs.MaxY(); })->MaxY();
+    double zMin = std::min_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MinZ() < rhs.MinZ(); })->MinZ();
+    double zMax = std::max_element(tpcs.begin(), tpcs.end(), [](auto &lhs, auto &rhs) { return lhs.MaxZ() < rhs.MaxZ(); })->MaxZ();
+    ActiveVolumes.emplace_back(xMin, xMax, yMin, yMax, zMin, zMax);
+  }
+  double fidInsetMinX = 0;
+  double fidInsetMaxX = 0;
+  double fidInsetMinY = 0;
+  double fidInsetMaxY = 0;
+  double fidInsetMinZ = 0;
+  double fidInsetMaxZ = 0;
+
+  if (fiducialVolumeInsets_.size() == 6) {
+    fidInsetMinX = fiducialVolumeInsets_[0];
+    fidInsetMaxX = fiducialVolumeInsets_[1];
+    fidInsetMinY = fiducialVolumeInsets_[2];
+    fidInsetMaxY = fiducialVolumeInsets_[3];
+    fidInsetMinZ = fiducialVolumeInsets_[4];
+    fidInsetMaxZ = fiducialVolumeInsets_[5];
+  }
+  else if (!fiducialVolumeInsets_.empty()) {
+    std::cout << "Error: fiducialVolumeInsets vector must have length of 6, not fiducializing" << std::endl;
+  }
+  for (const geo::BoxBoundedGeo &AV: ActiveVolumes) {
+    fiducialVolumes.emplace_back(AV.MinX() + fidInsetMinX, AV.MaxX() - fidInsetMaxX,
+                                 AV.MinY() + fidInsetMinY, AV.MaxY() - fidInsetMaxY,
+                                 AV.MinZ() + fidInsetMinZ, AV.MaxZ() - fidInsetMaxZ);
+  }
+  return fiducialVolumes;
+}
+
+std::vector<geo::BoxBoundedGeo> TrajectoryMCSFitter::setExcludeVolumes() const {
+  std::vector<geo::BoxBoundedGeo> excludeVolumes;
+  if (excludeVolumes_.size()%6 != 0) {
+    std::cout << "Error: excludeVolumes vector must have length multiple of 6, not excluding any regions" << std::endl;
+    excludeVolumes.emplace_back(-9999.0, -9999.0, -9999.0, -9999.0, -9999.0, -9999.0);
+    return excludeVolumes;
+  }
+  for (unsigned int i=0; i<excludeVolumes_.size()/6; i++) {
+    excludeVolumes.emplace_back(excludeVolumes_[6*i+0], excludeVolumes_[6*i+1], excludeVolumes_[6*i+2], excludeVolumes_[6*i+3], excludeVolumes_[6*i+4], excludeVolumes_[6*i+5]);
+  }
+  return excludeVolumes;
+}
+
+bool TrajectoryMCSFitter::isInVolume(const std::vector<geo::BoxBoundedGeo> &volumes, const geo::Point_t &point) const {
+  for (const geo::BoxBoundedGeo &volume: volumes) {
+    if (point.X()>=volume.MinX() && point.X()<=volume.MaxX() &&
+        point.Y()>=volume.MinY() && point.Y()<=volume.MaxY() &&
+        point.Z()>=volume.MinZ() && point.Z()<=volume.MaxZ()) {
+      return true;
+    }
+  }
+  return false;
+}

sbncode/LArRecoProducer/LArReco/TrajectoryMCSFitter.h

@@ -3,11 +3,15 @@
 
 #include "fhiclcpp/ParameterSet.h"
 #include "fhiclcpp/types/Atom.h"
+#include "fhiclcpp/types/Sequence.h"
 #include "fhiclcpp/types/Table.h"
 #include "canvas/Persistency/Common/Ptr.h"
 #include "lardataobj/RecoBase/MCSFitResult.h"
 #include "lardataobj/RecoBase/Track.h"
 #include "lardata/RecoObjects/TrackState.h"
+#include "larcore/Geometry/Geometry.h"
+#include "larcore/CoreUtils/ServiceUtil.h"
+#include "larcorealg/Geometry/GeometryCore.h"
 
 namespace trkf::sbn {
   /**
@@ -87,10 +91,18 @@ namespace trkf::sbn {
         Comment("Angular resolution parameter used in modified Highland formula. Unit is mrad."),
         3.0
       };
+      fhicl::Sequence<double> fiducialVolumeInsets {
+        Name("fiducialVolumeInsets"),
+        Comment("insets from the TPC boundaries to exclude from the fit")
+      };
+      fhicl::Sequence<double> excludeVolumes {
+        Name("excludeVolumes"),
+        Comment("regions to exclude")
+      };
     };
     using Parameters = fhicl::Table<Config>;
     //
-    TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol){
+    TrajectoryMCSFitter(int pIdHyp, int minNSegs, double segLen, int minHitsPerSegment, int nElossSteps, int eLossMode, double pMin, double pMax, double pStep, double angResol, std::vector<double>fiducialVolumeInsets, std::vector<double> excludeVolumes){
       pIdHyp_ = pIdHyp;
       minNSegs_ = minNSegs;
       segLen_ = segLen;
@@ -101,9 +113,11 @@ namespace trkf::sbn {
       pMax_ = pMax;
       pStep_ = pStep;
       angResol_ = angResol;
+      fiducialVolumeInsets_ = fiducialVolumeInsets;
+      excludeVolumes_ = excludeVolumes;
     }
     explicit TrajectoryMCSFitter(const Parameters & p)
-      : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol()) {}
+      : TrajectoryMCSFitter(p().pIdHypothesis(),p().minNumSegments(),p().segmentLength(),p().minHitsPerSegment(),p().nElossSteps(),p().eLossMode(),p().pMin(),p().pMax(),p().pStep(),p().angResol(),p().fiducialVolumeInsets(),p().excludeVolumes()) {}
     //
     recob::MCSFitResult fitMcs(const recob::TrackTrajectory& traj, bool momDepConst = true) const { return fitMcs(traj,pIdHyp_,momDepConst); }
     recob::MCSFitResult fitMcs(const recob::Track& track,          bool momDepConst = true) const { return fitMcs(track,pIdHyp_,momDepConst); }
@@ -117,7 +131,7 @@ namespace trkf::sbn {
       return fitMcs(tt,pid,momDepConst);
     }
     //
-    void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector<size_t>& breakpoints, std::vector<float>& segradlengths, std::vector<float>& cumseglens) const;
+    void breakTrajInSegments(const recob::TrackTrajectory& traj, std::vector<size_t>& breakpoints, std::vector<float>& segradlengths, std::vector<float>& cumseglens, std::vector<bool>& breakpointsgood) const;
     void linearRegression(const recob::TrackTrajectory& traj, const size_t firstPoint, const size_t lastPoint, recob::tracking::Vector_t& pcdir) const;
     double mcsLikelihood(double p, double theta0x, std::vector<float>& dthetaij, std::vector<float>& seg_nradl, std::vector<float>& cumLen, bool fwd, bool momDepConst, int pid) const;
     //
@@ -146,6 +160,9 @@ namespace trkf::sbn {
     double energyLossLandau(const double mass2,const double E2, const double x) const;
     //
     double GetE(const double initial_E, const double length_travelled, const double mass) const;
+    std::vector<geo::BoxBoundedGeo> setFiducialVolumes() const;
+    std::vector<geo::BoxBoundedGeo> setExcludeVolumes() const;
+    bool isInVolume(const std::vector<geo::BoxBoundedGeo> &volumes, const geo::Point_t &point) const;
     //
   private:
     int    pIdHyp_;
@@ -158,6 +175,8 @@ namespace trkf::sbn {
     double pMax_;
     double pStep_;
     double angResol_;
+    std::vector<double> fiducialVolumeInsets_;
+    std::vector<double> excludeVolumes_;
   };
 }
 

sbncode/LArRecoProducer/mcsproducer.fcl

@@ -1,7 +1,10 @@
 BEGIN_PROLOG
 mcs_sbn: {
   module_type: MCSFitAllPID
-  MCS: {}
+  MCS: {
+    fiducialVolumeInsets: [] 
+    excludeVolumes: [] 
+  }
   TrackLabel: pandoraTrack
 }
 

sbncode/LArRecoProducer/mcsproducer_icarus.fcl

@@ -0,0 +1,13 @@
+BEGIN_PROLOG
+mcs_sbn: {
+  module_type: MCSFitAllPID
+  MCS: {
+    eLossMode: 2
+    angResol: 10.0
+    fiducialVolumeInsets: [10.0, 10.0, 10.0, 10.0, 10.0, 10.0] 
+    excludeVolumes: [190.0, 240.0, -9999.0, 9999.0, -9999.0, 9999.0] 
+  }
+  TrackLabel: pandoraTrack
+}
+
+END_PROLOG