Feature/jlarkin numu selection

sbnana/SBNAna/Analysis/202208Scanning/selected_events.C

@@ -0,0 +1,38 @@
+#include "sbnana/CAFAna/Core/Binning.h"
+#include "sbnana/CAFAna/Core/Cut.h"
+#include "sbnana/CAFAna/Core/Spectrum.h"
+#include "sbnana/CAFAna/Core/SpectrumLoader.h"
+#include "sbnana/CAFAna/Core/Var.h"
+
+#include "sbnana/SBNAna/Cuts/NumuCutsIcarus202208.h"
+#include "sbnana/SBNAna/Vars/NumuVarsIcarus202208.h"
+
+#include "sbnanaobj/StandardRecord/Proxy/SRProxy.h"
+
+#include <fstream>
+
+using namespace ana;
+
+void selected_events()
+{
+  SpectrumLoader loader("icarus_BNB_Nu_Cosmics_v09_37_02_04_flatcaf");
+
+  std::ofstream out("selected_events.txt");
+
+  const SpillVar dummy_var([&out](const caf::SRSpillProxy* spill){
+    for(const auto& slc: spill->slc) {
+      if(kIcarus202208QELikeContainedMuon(&slc)) {
+        out << "Run: " << spill->hdr.run << "\tSubrun: " << spill->hdr.subrun
+            << "\tEvent: " << spill->hdr.evt << "\tVertex: (" 
+            << slc.vertex.x << ", " << slc.vertex.y << ", " << slc.vertex.z << ")\n";
+        return 1;
+      }
+    } 
+    return 0;
+  });
+
+  const Binning dummy_bins = Binning::Simple(2, 0, 2);
+  Spectrum dummy_spec("Dummy Label", dummy_bins, loader, dummy_var, kNoSpillCut);
+
+  loader.Go();
+}

sbnana/SBNAna/Cuts/NumuCutsIcarus202208.cxx

@@ -0,0 +1,51 @@
+#include "sbnana/CAFAna/Core/Cut.h"
+#include "sbnana/SBNAna/Cuts/NumuCutsIcarus202208.h"
+#include "sbnana/SBNAna/Vars/NumuVarsIcarus202208.h"
+
+#include "sbnanaobj/StandardRecord/Proxy/FwdDeclare.h"
+#include "sbnanaobj/StandardRecord/Proxy/SRProxy.h"
+
+namespace ana {
+
+static bool contained(const caf::SRTrackProxy& trk)
+{
+  return ((trk.end.x < -71.1 - 5 && trk.end.x > -369.33 + 5)
+             || (trk.end.x > 71.1 + 5 && trk.end.x < 369.33 - 5))
+           && trk.end.y > -181.7 + 5 && trk.end.y < 134.8 - 5
+           && trk.end.z > -895.95 + 5 && trk.end.z < 895.95 - 5;
+}
+
+const Cut kIcarus202208FMTimeCut = kFMTimeVar > 0 && kFMTimeVar < 1.8;
+const Cut kIcarus202208FMScoreCut = kFMScoreVar < 9;
+const Cut kIcarus202208LongTrackDirCut = kCRLongestTrackDirY > -0.91;
+const Cut kIcarus202208FoundMuon = kIcarus202208MuonIdx >= 0;
+const Cut kIcarus202208RecoFiducial([](const caf::SRSliceProxy* slc) {
+      return ( !isnan(slc->vertex.x) &&
+	       ( ( slc->vertex.x < -71.1 - 25 && slc->vertex.x > -369.33 + 25 ) ||
+		 ( slc->vertex.x > 71.1 + 25 && slc->vertex.x < 369.33 - 25 ) ) &&
+	       !isnan(slc->vertex.y) &&
+	       ( slc->vertex.y > -181.7 + 25 && slc->vertex.y < 134.8 - 25 ) &&
+	       !isnan(slc->vertex.z) &&
+	       ( slc->vertex.z > -895.95 + 30 && slc->vertex.z < 895.95 - 50 ) );
+    });
+
+const Cut kIcarus202208NumuSelection = kIcarus202208RecoFiducial && kIcarus202208FMScoreCut && kIcarus202208FMTimeCut && kIcarus202208LongTrackDirCut && kIcarus202208FoundMuon;
+
+const Cut kIcarus202208NoPion = kIcarus202208NumPions == 0;
+const Cut kIcarus202208ContainedHadrons([](const caf::SRSliceProxy* slc){
+  auto idx = kIcarus202208MuonIdx(slc);
+  int muID = -1;
+  if (idx >= 0) muID = slc->reco.trk.at(idx).pfp.id;
+  for(auto& trk: slc->reco.trk) {
+    if(trk.pfp.id != muID && trk.pfp.parent_is_primary)
+      if(!contained(trk)) return false;
+  }
+  return true;
+});
+const Cut kIcarus202208ContainedMuon([](const caf::SRSliceProxy* slc){
+  return kIcarus202208FoundMuon(slc) && contained(slc->reco.trk.at(kIcarus202208MuonIdx(slc)));
+});
+const Cut kIcarus202208ContainedMuonAndHadrons = kIcarus202208ContainedMuon && kIcarus202208ContainedHadrons;
+const Cut kIcarus202208QELike = kIcarus202208NumuSelection && kIcarus202208NoPion && kIcarus202208ContainedHadrons;
+const Cut kIcarus202208QELikeContainedMuon = kIcarus202208QELike && kIcarus202208ContainedMuon;
+}

sbnana/SBNAna/Cuts/NumuCutsIcarus202208.h

@@ -0,0 +1,21 @@
+#pragma once
+
+#include "sbnana/CAFAna/Core/Cut.h"
+
+namespace ana {
+
+extern const Cut kIcarus202208FMTimeCut;
+extern const Cut kIcarus202208FMScoreCut;
+extern const Cut kIcarus202208LongTrackDirCut;
+extern const Cut kIcarus202208FoundMuon;
+extern const Cut kIcarus202208RecoFiducial;
+
+extern const Cut kIcarus202208NumuSelection;
+
+extern const Cut kIcarus202208NoPion;
+extern const Cut kIcarus202208ContainedHadrons;
+extern const Cut kIcarus202208ContainedMuon;
+extern const Cut kIcarus202208ContainedMuonAndHadrons;
+extern const Cut kIcarus202208QELike;
+extern const Cut kIcarus202208QELikeContainedMuon;
+}

sbnana/SBNAna/Vars/NumuVarsIcarus202208.cxx

@@ -0,0 +1,65 @@
+#include "sbnana/SBNAna/Vars/NumuVarsIcarus202208.h"
+#include "sbnana/CAFAna/Core/Var.h"
+
+#include "sbnanaobj/StandardRecord/Proxy/FwdDeclare.h"
+#include "sbnanaobj/StandardRecord/Proxy/SRProxy.h"
+
+namespace ana {
+
+// More general vars, i.e., not just numu vars.
+// Probably should be moved to another file at some point.
+const Var kCRLongestTrackDirY = SIMPLEVAR(nuid.crlongtrkdiry);
+const Var kFMTimeVar          = SIMPLEVAR(fmatch.time);
+const Var kFMScoreVar         = SIMPLEVAR(fmatch.score);
+
+// Mostly just copied from NumuVarsIcarus202106.cxx
+const Var kIcarus202208MuonIdx([](const caf::SRSliceProxy* slc) -> int {
+      // The (dis)qualification of a slice is based upon the track level information.
+      float Longest(0);
+      int PTrackInd(-1);
+      for (std::size_t i(0); i < slc->reco.trk.size(); ++i) {
+        auto const& trk = slc->reco.trk.at(i);
+        if(trk.bestplane == -1) continue;
+        const float Atslc = std::hypot(slc->vertex.x - trk.start.x,
+                                       slc->vertex.y - trk.start.y,
+                                       slc->vertex.z - trk.start.z);
+        const bool AtSlice = ( Atslc < 10.0 && trk.pfp.parent_is_primary);
+
+        const float Chi2Proton = trk.chi2pid[trk.bestplane].chi2_proton;
+        const float Chi2Muon = trk.chi2pid[trk.bestplane].chi2_muon;
+
+        const bool Contained = ( !isnan(trk.end.x) &&
+        (( trk.end.x < -71.1 - 5 && trk.end.x > -369.33 + 5 ) ||
+        ( trk.end.x < 71.1 + 5 && trk.end.x > +369.33 - 5 )) &&
+        !isnan(trk.end.y) &&
+        ( trk.end.y > -181.7 + 5 && trk.end.y < 134.8 - 5 ) &&
+        !isnan(trk.end.z) &&
+        ( trk.end.z > -895.95 + 5 && trk.end.z < 895.95 - 5 ) );
+        const bool MaybeMuonExiting = ( !Contained && trk.len > 100);
+        const bool MaybeMuonContained = ( Contained && Chi2Proton > 60 && Chi2Muon < 30 && trk.len > 50 );
+        if ( AtSlice && ( MaybeMuonExiting || MaybeMuonContained ) && trk.len > Longest ) {
+          Longest = trk.len;
+          PTrackInd = i;
+        }
+      }
+      return PTrackInd;
+});
+
+static bool proton_cut(const caf::SRTrackProxy& trk)
+{
+  return trk.chi2pid[trk.bestplane].chi2_proton < 100;
+}
+
+const Var kIcarus202208NumPions([](const caf::SRSliceProxy* slc)
+{
+  int count = 0;
+  auto idx = kIcarus202208MuonIdx(slc);
+  int muID = -1;
+  if (idx >= 0) muID = slc->reco.trk.at(idx).pfp.id;
+  for(auto& trk: slc->reco.trk) {
+    if(trk.pfp.id != muID && !proton_cut(trk) && trk.pfp.parent_is_primary)
+      ++count;
+  }
+  return count;
+});
+}

sbnana/SBNAna/Vars/NumuVarsIcarus202208.h

@@ -0,0 +1,15 @@
+#pragma once
+
+#include "sbnana/CAFAna/Core/Var.h"
+
+namespace ana {
+
+extern const Var kCRLongestTrackDirY;
+extern const Var kFMTimeVar;
+extern const Var kFMScoreVar;
+
+extern const Var kIcarus202208MuonIdx;
+extern const Var kIcarus202208NumPions;
+}
+
+