diff --git a/fcl/g4/larg4_icarus_cosmics_sce_filterdirt.fcl b/fcl/g4/larg4_icarus_cosmics_sce_filterdirt.fcl new file mode 100644 index 000000000..f4c0d3a52 --- /dev/null +++ b/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" diff --git a/fcl/gen/genie/simulation_genie_icarus_bnb_volDetEnclosure_fiducial.fcl b/fcl/gen/genie/simulation_genie_icarus_bnb_volDetEnclosure_fiducial.fcl new file mode 100644 index 000000000..09612bf48 --- /dev/null +++ b/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" diff --git a/fcl/gen/genie/simulation_genie_icarus_bnb_volWorld_DetEnclosurerockbox.fcl b/fcl/gen/genie/simulation_genie_icarus_bnb_volWorld_DetEnclosurerockbox.fcl new file mode 100644 index 000000000..aa800408b --- /dev/null +++ b/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"] diff --git a/fcl/gen/genie/simulation_genie_icarus_bnb_volWorld_rockbox.fcl b/fcl/gen/genie/simulation_genie_icarus_bnb_volWorld_rockbox.fcl new file mode 100644 index 000000000..f8a411f55 --- /dev/null +++ b/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/" diff --git a/fcl/overlays/dirt_gen_overlay_siminfomixer.fcl b/fcl/overlays/dirt_gen_overlay_siminfomixer.fcl new file mode 100644 index 000000000..ebb55b3dd --- /dev/null +++ b/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"] diff --git a/icaruscode/Filters/FilterMCTruthVolume_module.cc b/icaruscode/Filters/FilterMCTruthVolume_module.cc new file mode 100644 index 000000000..8aac0b1fa --- /dev/null +++ b/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 +#include +#include + +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 fVolumeLow; + std::array fVolumeHigh; + + bool inVolume(double x, double y, double z) const; +}; + +// --------------------------------------------------------------------------- +FilterMCTruthVolume::FilterMCTruthVolume(fhicl::ParameterSet const& pset) + : EDFilter{pset} + , fMCTruthLabel{pset.get("MCTruthLabel", "generator")} + , fPOTLabel{pset.get("POTLabel", "generator")} +{ + auto lo = pset.get>("VolumeLow"); + auto hi = pset.get>("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>(); + produces>(); + produces>(); + produces>(); + produces>(); + produces(); +} + +// --------------------------------------------------------------------------- +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>(); + auto outGTruth = std::make_unique>(); + auto outMCFlux = std::make_unique>(); + auto assnsTG = std::make_unique>(); + auto assnsTF = std::make_unique>(); + + auto const& mctruthHandle = evt.getValidHandle>(fMCTruthLabel); + + art::FindOneP findGTruth(mctruthHandle, evt, fMCTruthLabel); + art::FindOneP findMCFlux(mctruthHandle, evt, fMCTruthLabel); + + art::PtrMaker makeTruthPtr(evt); + art::PtrMaker makeGTruthPtr(evt); + art::PtrMaker 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(fPOTLabel); + sr.put(std::make_unique(*potHandle), art::subRunFragment()); + return true; +} + +DEFINE_ART_MODULE(FilterMCTruthVolume) diff --git a/icaruscode/Filters/SimEnergyDepFakeTriggerFilterICARUS_module.cc b/icaruscode/Filters/SimEnergyDepFakeTriggerFilterICARUS_module.cc new file mode 100644 index 000000000..b986d7a40 --- /dev/null +++ b/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("BeamTimeMin")) + , fBeamTimeMax(pset.get("BeamTimeMax")) + , fEnergyDeposit(pset.get("EnergyDeposit")) + , fMaxEnergyDeposit(pset.get("MaxEnergyDeposit", std::numeric_limits::max())) + , fSimEnergyDepModuleName(pset.get("SimEnergyDepModuleName")) +{ +} + +bool SimEnergyDepFakeTriggerFilterICARUS::filter(art::Event& e) +{ + const art::ValidHandle>& + energyDeps(e.getValidHandle>(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) + +} + diff --git a/icaruscode/Filters/filterMCTruthVolume.fcl b/icaruscode/Filters/filterMCTruthVolume.fcl new file mode 100644 index 000000000..32978a7cd --- /dev/null +++ b/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 diff --git a/icaruscode/Filters/simenergydep_faketrigger_icarus.fcl b/icaruscode/Filters/simenergydep_faketrigger_icarus.fcl new file mode 100644 index 000000000..535c5ba2f --- /dev/null +++ b/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 diff --git a/icaruscode/Overlays/SubRunPOTInEvent_module.cc b/icaruscode/Overlays/SubRunPOTInEvent_module.cc new file mode 100644 index 000000000..60a97a2db --- /dev/null +++ b/icaruscode/Overlays/SubRunPOTInEvent_module.cc @@ -0,0 +1,90 @@ +//////////////////////////////////////////////////////////////////////// +// Class: SubRunPOTInEvent +// Plugin Type: producer (art v3_00_00) +// File: SubRunPOTInEvent_module.cc +// +// Generated at Mon Jan 21 09:09:38 2019 by Wesley Ketchum using cetskelgen +// from cetlib version v3_04_00. +// +// This module directly copies the subrun POT info into every event. +// **IT IS NOT A POT PER EVENT IT IS NOT A POT PER EVENT** +// **IT IS NOT A POT PER EVENT IT IS NOT A POT PER EVENT** +// **IT IS NOT A POT PER EVENT IT IS NOT A POT PER EVENT** +// **IT IS NOT A POT PER EVENT IT IS NOT A POT PER EVENT** +// **IT IS NOT A POT PER EVENT IT IS NOT A POT PER EVENT** +// **IT IS NOT A POT PER EVENT IT IS NOT A POT PER EVENT** +// +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include + +#include "larcoreobj/SummaryData/POTSummary.h" + +namespace mix { + class SubRunPOTInEvent; +} + + +class mix::SubRunPOTInEvent : public art::EDProducer { +public: + explicit SubRunPOTInEvent(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + SubRunPOTInEvent(SubRunPOTInEvent const&) = delete; + SubRunPOTInEvent(SubRunPOTInEvent&&) = delete; + SubRunPOTInEvent& operator=(SubRunPOTInEvent const&) = delete; + SubRunPOTInEvent& operator=(SubRunPOTInEvent&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + + // Selected optional functions. + void beginSubRun(art::SubRun& sr) override; + +private: + + art::InputTag fInputTag; + sumdata::POTSummary fSubRunPOT; + +}; + + +mix::SubRunPOTInEvent::SubRunPOTInEvent(fhicl::ParameterSet const& p) + : EDProducer{p} // , + // More initializers here. +{ + fInputTag = p.get("InputTag"); + + produces("SubRunPOT"); +} + +void mix::SubRunPOTInEvent::produce(art::Event& e) +{ + std::unique_ptr srpot_ptr(new sumdata::POTSummary(fSubRunPOT)); + e.put(std::move(srpot_ptr),"SubRunPOT"); +} + +void mix::SubRunPOTInEvent::beginSubRun(art::SubRun& sr) +{ + art::Handle srpot_handle; + sr.getByLabel(fInputTag,srpot_handle); + if(srpot_handle.isValid()) fSubRunPOT = *srpot_handle; + else std::cout << "BAD SUBRUNPOT\n"; + + std::cout << "GOT POT: " << fSubRunPOT.totgoodpot << std::endl; +} + +DEFINE_ART_MODULE(mix::SubRunPOTInEvent) +