New fcl configurations and modules to enable Dirt Overlay configuration for ICARUS Runs 2+4 MC.

fcl/g4/larg4_icarus_cosmics_sce_filterdirt.fcl

@@ -0,0 +1,11 @@
+#include "simenergydep_faketrigger_icarus.fcl"
+#include "icarus_siminfomixer.fcl"
+#include "larg4_icarus_cosmics.fcl"
+
+physics.filters.dirtfilter: @local::icarus_simenergydepfaketriggerfilter
+physics.producers.potinevent: @local::icarus_subrunpotinevent # so the sim info mixer can grab it
+
+outputs.rootoutput.SelectEvents: [simulate]
+physics.simulate: [@sequence::physics.simulate, potinevent, dirtfilter]
+
+#include "enable_spacecharge_icarus.fcl"

fcl/gen/genie/simulation_genie_icarus_bnb_volDetEnclosure_fiducial.fcl

@@ -0,0 +1,3 @@
+#include "simulation_genie_icarus_bnb_volDetEnclosure.fcl"
+
+physics.producers.generator.FiducialCut: "mbox: -378.49,-191.86,-904.950652270838,378.49,144.96,904.950652270838"

fcl/gen/genie/simulation_genie_icarus_bnb_volWorld_DetEnclosurerockbox.fcl

@@ -0,0 +1,21 @@
+#include "filterMCTruthVolume.fcl"
+#include "simulation_genie_icarus_bnb_volDetEnclosure.fcl"
+
+# change "generator" to "mcgen"-- we'll have a filter module to replicate generator
+physics.producers.mcgen: @local::physics.producers.generator
+physics.producers.mcgen.TopVolume:  "volWorld"
+physics.producers.mcgen.FiducialCut: "rockbox:(-378.49,-191.86,-904.950652270838)(+378.49,+144.96,+904.950652270838),1,800,0.00425,1.3,1"
+physics.producers.mcgen.FluxFiles: ["gsimple_bnb_neutrino_icarus_dirt_*.root"]
+physics.producers.mcgen.FluxSearchPaths: "/cvmfs/sbn.osgstorage.org/pnfs/fnal.gov/usr/sbn/persistent/stash/physics-gputnam/icarus-bnb-dirt/"
+physics.producers.generator: @erase
+
+# Reject neutrino interactions inside volDetEnclosure
+physics.filters.generator: @local::filtermctruthvolume
+
+physics.simulate: [rns, mcgen, generator, beamgate]
+
+outputs.rootoutput.outputCommands: [
+  "keep *_*_*_*",
+  "drop *_mcgen_*_*"
+]
+outputs.rootoutput.SelectEvents: ["simulate"]

fcl/gen/genie/simulation_genie_icarus_bnb_volWorld_rockbox.fcl

@@ -0,0 +1,6 @@
+#include "simulation_genie_icarus_bnb_volDetEnclosure.fcl"
+
+physics.producers.generator.TopVolume:  "volWorld"
+physics.producers.generator.FiducialCut: "rockbox:(-378.49,-191.86,-904.950652270838)(+378.49,+144.96,+904.950652270838),1,800,0.00425,1.3,1"
+physics.producers.generator.FluxFiles: ["gsimple_bnb_neutrino_icarus_dirt_*.root"]
+physics.producers.generator.FluxSearchPaths: "/cvmfs/sbn.osgstorage.org/pnfs/fnal.gov/usr/sbn/persistent/stash/physics-gputnam/icarus-bnb-dirt/"

fcl/overlays/dirt_gen_overlay_siminfomixer.fcl

@@ -0,0 +1,111 @@
+#include "services_common_icarus.fcl"
+#include "icarus_siminfomixer.fcl"
+
+services:
+{
+  @table::icarus_common_services
+
+  scheduler:    { defaultExceptions: false }    # Make all uncaught exceptions fatal.
+  # Load the service that manages root files for histograms.
+  TFileService: { fileName: "siminfomixer.root" }
+
+  RandomNumberGenerator: {} #ART native random number generator
+}
+
+process_name : SimInfoMixer     #The process name must NOT contain any underscores
+
+source:
+{
+  module_type: RootInput
+  saveMemoryObjectThreshold: 0
+  maxEvents: -1
+}
+
+outputs: {
+  out: { module_type: RootOutput
+  	 fileName: "%ifb_%tc_mixed.root"
+	 compressionLevel: 1
+	 dataTier: "generated"
+         SelectEvents: ["mixer_path"]
+       }
+}
+
+physics: {
+
+  producers : {
+  }
+
+  analyzers: {
+  }
+
+  filters : {
+    generator: @local::icarus_siminfomixer
+    beamgate: @local::icarus_siminfomixer
+    largeant: @local::icarus_siminfomixer
+    ionization: @local::icarus_siminfomixer
+    simplemerge: @local::icarus_siminfomixer
+    genericcrt: @local::icarus_siminfomixer
+    sedlite: @local::icarus_siminfomixer
+    simdrift: @local::icarus_siminfomixer
+    pdfastsim: @local::icarus_siminfomixer
+  }
+
+  mixer_path : [ generator, beamgate, largeant, ionization, simplemerge, genericcrt, sedlite, simdrift, pdfastsim]
+  trigger_paths : [ mixer_path ]
+
+  output : [ out ]
+  end_paths: [ output ]
+
+}
+
+physics.filters.generator.MCTruthInputModuleLabels: ["generator::GenGenie"]
+physics.filters.generator.GTruthInputModuleLabels: ["generator::GenGenie"]
+physics.filters.generator.MCTruthGTruthAssnsInputModuleLabels: ["generator::GenGenie"]
+physics.filters.generator.MCFluxInputModuleLabels: ["generator::GenGenie"]
+physics.filters.generator.MCTruthMCFluxAssnsInputModuleLabels: ["generator::GenGenie"]
+physics.filters.generator.FillPOTInfo: true
+physics.filters.generator.POTSummaryTag: "potinevent:SubRunPOT:G4"
+
+physics.filters.beamgate.BeamGateInputModuleLabels: ["beamgate::GenGenie"]
+
+physics.filters.largeant.SimEnergyDepositInputModuleLabels: [
+  "largeant:LArG4DetectorServicevolTPCPlaneV:G4",
+  "largeant:LArG4DetectorServicevolTPC0:G4",
+  "largeant:LArG4DetectorServicevolTPCPlaneU:G4",
+  "largeant:LArG4DetectorServicevolTPCPlaneY:G4",
+  "largeant:LArG4DetectorServicevolTPCActive:G4"
+]
+physics.filters.largeant.AuxDetHitInputModuleLabels: [
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveCERNbot:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut309:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOS:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut256:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveDC:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut508:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut497:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut325:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut485:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveCERNtop:G4",
+  "largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut400:G4"
+]
+physics.filters.largeant.MCParticleInputModuleLabels: [
+  "largeant::G4",
+  "largeant:droppedMCParticles:G4"
+]
+physics.filters.largeant.MCTruthInputModuleLabels: ["generator::GenGenie"]
+physics.filters.largeant.MCTruthMCParticleAssnsInputModuleLabels: ["largeant::G4"]
+
+physics.filters.ionization.SimEnergyDepositInputModuleLabels: [
+  "ionization::G4",
+  "ionization:priorSCE:G4"
+]
+
+physics.filters.simplemerge.MCParticleInputModuleLabels: ["simplemerge::G4"]
+
+physics.filters.genericcrt.AuxDetSimChannelInputModuleLabels: ["genericcrt::G4"]
+
+physics.filters.sedlite.SimEnergyDepositLiteInputModuleLabels: ["sedlite::G4"]
+
+physics.filters.simdrift.SimChannelInputModuleLabels: ["simdrift::G4"]
+
+physics.filters.pdfastsim.SimPhotonsInputModuleLabels: ["pdfastsim::G4"]

icaruscode/Filters/FilterMCTruthVolume_module.cc

@@ -0,0 +1,145 @@
+////////////////////////////////////////////////////////////////////////
+// FilterMCTruthVolume_module.cc
+//
+// art::EDFilter that selects events containing at least one MCTruth
+// neutrino interaction with a vertex __outside__ a user-specified box.
+// Produces filtered MCTruth, GTruth, and MCFlux collections with
+// their associations. Passes through POT summary at the SubRun level.
+////////////////////////////////////////////////////////////////////////
+
+#include "art/Framework/Core/EDFilter.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Principal/SubRun.h"
+#include "art/Persistency/Common/PtrMaker.h"
+#include "canvas/Persistency/Common/Assns.h"
+#include "canvas/Persistency/Common/FindOneP.h"
+#include "canvas/Utilities/InputTag.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "larcoreobj/SummaryData/POTSummary.h"
+#include "nusimdata/SimulationBase/MCTruth.h"
+#include "nusimdata/SimulationBase/GTruth.h"
+#include "nusimdata/SimulationBase/MCFlux.h"
+
+#include <array>
+#include <memory>
+#include <vector>
+
+class FilterMCTruthVolume : public art::EDFilter {
+public:
+  explicit FilterMCTruthVolume(fhicl::ParameterSet const& pset);
+
+  bool filter(art::Event& evt) override;
+  bool endSubRun(art::SubRun& sr) override;
+
+private:
+  art::InputTag fMCTruthLabel;
+  art::InputTag fPOTLabel;
+  std::array<double, 3> fVolumeLow;
+  std::array<double, 3> fVolumeHigh;
+
+  bool inVolume(double x, double y, double z) const;
+};
+
+// ---------------------------------------------------------------------------
+FilterMCTruthVolume::FilterMCTruthVolume(fhicl::ParameterSet const& pset)
+  : EDFilter{pset}
+  , fMCTruthLabel{pset.get<art::InputTag>("MCTruthLabel", "generator")}
+  , fPOTLabel{pset.get<art::InputTag>("POTLabel", "generator")}
+{
+  auto lo = pset.get<std::vector<double>>("VolumeLow");
+  auto hi = pset.get<std::vector<double>>("VolumeHigh");
+  if (lo.size() != 3 || hi.size() != 3)
+    throw art::Exception(art::errors::Configuration)
+      << "VolumeLow and VolumeHigh must each have exactly 3 elements (x, y, z).";
+  std::copy(lo.begin(), lo.end(), fVolumeLow.begin());
+  std::copy(hi.begin(), hi.end(), fVolumeHigh.begin());
+
+  produces<std::vector<simb::MCTruth>>();
+  produces<std::vector<simb::GTruth>>();
+  produces<std::vector<simb::MCFlux>>();
+  produces<art::Assns<simb::MCTruth, simb::GTruth>>();
+  produces<art::Assns<simb::MCTruth, simb::MCFlux>>();
+  produces<sumdata::POTSummary, art::InSubRun>();
+}
+
+// ---------------------------------------------------------------------------
+bool FilterMCTruthVolume::inVolume(double x, double y, double z) const
+{
+  return x >= fVolumeLow[0] && x <= fVolumeHigh[0] &&
+         y >= fVolumeLow[1] && y <= fVolumeHigh[1] &&
+         z >= fVolumeLow[2] && z <= fVolumeHigh[2];
+}
+
+// ---------------------------------------------------------------------------
+bool FilterMCTruthVolume::filter(art::Event& evt)
+{
+  auto outMCTruth = std::make_unique<std::vector<simb::MCTruth>>();
+  auto outGTruth  = std::make_unique<std::vector<simb::GTruth>>();
+  auto outMCFlux  = std::make_unique<std::vector<simb::MCFlux>>();
+  auto assnsTG    = std::make_unique<art::Assns<simb::MCTruth, simb::GTruth>>();
+  auto assnsTF    = std::make_unique<art::Assns<simb::MCTruth, simb::MCFlux>>();
+
+  auto const& mctruthHandle = evt.getValidHandle<std::vector<simb::MCTruth>>(fMCTruthLabel);
+
+  art::FindOneP<simb::GTruth> findGTruth(mctruthHandle, evt, fMCTruthLabel);
+  art::FindOneP<simb::MCFlux> findMCFlux(mctruthHandle, evt, fMCTruthLabel);
+
+  art::PtrMaker<simb::MCTruth> makeTruthPtr(evt);
+  art::PtrMaker<simb::GTruth>  makeGTruthPtr(evt);
+  art::PtrMaker<simb::MCFlux>  makeMCFluxPtr(evt);
+
+  for (size_t i = 0; i < mctruthHandle->size(); ++i) {
+    auto const& mct = mctruthHandle->at(i);
+
+    if (!mct.NeutrinoSet()) continue;
+
+    double vx = mct.GetNeutrino().Nu().Vx();
+    double vy = mct.GetNeutrino().Nu().Vy();
+    double vz = mct.GetNeutrino().Nu().Vz();
+
+    if (inVolume(vx, vy, vz)) continue;
+
+    outMCTruth->push_back(mct);
+    size_t idx = outMCTruth->size() - 1;
+
+    if (findGTruth.isValid()) {
+      auto const& gtp = findGTruth.at(i);
+      if (gtp.isNonnull()) {
+        outGTruth->push_back(*gtp);
+        assnsTG->addSingle(makeTruthPtr(idx), makeGTruthPtr(outGTruth->size() - 1));
+      }
+    }
+
+    if (findMCFlux.isValid()) {
+      auto const& mfp = findMCFlux.at(i);
+      if (mfp.isNonnull()) {
+        outMCFlux->push_back(*mfp);
+        assnsTF->addSingle(makeTruthPtr(idx), makeMCFluxPtr(outMCFlux->size() - 1));
+      }
+    }
+  }
+
+  bool pass = !outMCTruth->empty();
+
+  evt.put(std::move(outMCTruth));
+  evt.put(std::move(outGTruth));
+  evt.put(std::move(outMCFlux));
+  evt.put(std::move(assnsTG));
+  evt.put(std::move(assnsTF));
+
+  return pass;
+}
+
+// ---------------------------------------------------------------------------
+bool FilterMCTruthVolume::endSubRun(art::SubRun& sr)
+{
+  auto const& potHandle = sr.getValidHandle<sumdata::POTSummary>(fPOTLabel);
+  sr.put(std::make_unique<sumdata::POTSummary>(*potHandle), art::subRunFragment());
+  return true;
+}
+
+DEFINE_ART_MODULE(FilterMCTruthVolume)

icaruscode/Filters/SimEnergyDepFakeTriggerFilterICARUS_module.cc

@@ -0,0 +1,58 @@
+#include "art/Framework/Core/EDFilter.h"
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Principal/Event.h"
+
+#include "lardataobj/Simulation/SimEnergyDeposit.h"
+
+namespace filt {
+
+class SimEnergyDepFakeTriggerFilterICARUS : public art::EDFilter {
+  public:
+  explicit SimEnergyDepFakeTriggerFilterICARUS(fhicl::ParameterSet const& pset);
+  virtual bool filter(art::Event& e) override;
+
+  private:
+  const double fBeamTimeMin;   // Minimum time of beam window [us]
+  const double fBeamTimeMax;   // Maximum time of beam window [us]
+  const double fEnergyDeposit; // Minimum energy deposit in TPC for trigger [MeV]
+  const double fMaxEnergyDeposit; // Maximum energy deposit in TPC for trigger [MeV]
+
+  const std::string fSimEnergyDepModuleName;
+};
+
+SimEnergyDepFakeTriggerFilterICARUS::SimEnergyDepFakeTriggerFilterICARUS(fhicl::ParameterSet const& pset)
+    : EDFilter(pset)
+    , fBeamTimeMin(pset.get<double>("BeamTimeMin"))
+    , fBeamTimeMax(pset.get<double>("BeamTimeMax"))
+    , fEnergyDeposit(pset.get<double>("EnergyDeposit"))
+    , fMaxEnergyDeposit(pset.get<double>("MaxEnergyDeposit", std::numeric_limits<double>::max()))
+    , fSimEnergyDepModuleName(pset.get<std::string>("SimEnergyDepModuleName"))
+{
+}
+
+bool SimEnergyDepFakeTriggerFilterICARUS::filter(art::Event& e)
+{
+  const art::ValidHandle<std::vector<sim::SimEnergyDeposit>>&
+      energyDeps(e.getValidHandle<std::vector<sim::SimEnergyDeposit>>(fSimEnergyDepModuleName));
+
+  double energy(0);
+
+  for (const sim::SimEnergyDeposit& energyDep : *energyDeps) {
+    // Check particle time is within the beam time
+    const double time(energyDep.Time() * 1e-3); // [ns] -> [us]
+    if (time < fBeamTimeMin || time > fBeamTimeMax)
+      continue;
+
+    // Add up the energy deposit inside the TPC
+    energy += energyDep.Energy(); // [MeV]
+  }
+
+  std::cout << "SAW E= " << energy << " inside window. Threshold is: " << fEnergyDeposit << std::endl;
+  // If the energy deposit within the beam time is greater than some limit then trigger the event
+  return (energy > fEnergyDeposit) & (energy < fMaxEnergyDeposit);
+}
+
+DEFINE_ART_MODULE(SimEnergyDepFakeTriggerFilterICARUS)
+
+}
+

icaruscode/Filters/filterMCTruthVolume.fcl

@@ -0,0 +1,11 @@
+BEGIN_PROLOG
+
+filtermctruthvolume: {
+  module_type: FilterMCTruthVolume
+  MCTruthLabel: "mcgen"
+  POTLabel: "mcgen"
+  VolumeLow:  [-571.55, -352.067, -1187.04]
+  VolumeHigh: [+571.55, +619.033, +1566.2]
+}
+
+END_PROLOG

icaruscode/Filters/simenergydep_faketrigger_icarus.fcl

@@ -0,0 +1,15 @@
+BEGIN_PROLOG
+
+icarus_simenergydepfaketriggerfilter:
+{
+  module_type:     "SimEnergyDepFakeTriggerFilterICARUS"
+
+  BeamTimeMin:    -0.2        # Minimum time of beam window [us]
+  BeamTimeMax:     1.9        # Maximum time of beam window [us]
+  EnergyDeposit:   5          # Minimum energy deposit in TPC for trigger [MeV]
+
+  # By default, take only the energy deposits within the TPC active volume
+  SimEnergyDepModuleName: "largeant:LArG4DetectorServicevolTPCActive" # Name of SimEnergyDeposit producer module
+}
+
+END_PROLOG