Skip to content

Commit 9ec6034

Browse files
authored
Merge pull request #894 from SBNSoftware/feature/gputnam-dirt-workflows
New fcl configurations and modules to enable Dirt Overlay configuration for ICARUS Runs 2+4 MC.
2 parents 47f5a15 + 3232ba3 commit 9ec6034

10 files changed

Lines changed: 471 additions & 0 deletions
Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
#include "simenergydep_faketrigger_icarus.fcl"
2+
#include "icarus_siminfomixer.fcl"
3+
#include "larg4_icarus_cosmics.fcl"
4+
5+
physics.filters.dirtfilter: @local::icarus_simenergydepfaketriggerfilter
6+
physics.producers.potinevent: @local::icarus_subrunpotinevent # so the sim info mixer can grab it
7+
8+
outputs.rootoutput.SelectEvents: [simulate]
9+
physics.simulate: [@sequence::physics.simulate, potinevent, dirtfilter]
10+
11+
#include "enable_spacecharge_icarus.fcl"
Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
#include "simulation_genie_icarus_bnb_volDetEnclosure.fcl"
2+
3+
physics.producers.generator.FiducialCut: "mbox: -378.49,-191.86,-904.950652270838,378.49,144.96,904.950652270838"
Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
#include "filterMCTruthVolume.fcl"
2+
#include "simulation_genie_icarus_bnb_volDetEnclosure.fcl"
3+
4+
# change "generator" to "mcgen"-- we'll have a filter module to replicate generator
5+
physics.producers.mcgen: @local::physics.producers.generator
6+
physics.producers.mcgen.TopVolume: "volWorld"
7+
physics.producers.mcgen.FiducialCut: "rockbox:(-378.49,-191.86,-904.950652270838)(+378.49,+144.96,+904.950652270838),1,800,0.00425,1.3,1"
8+
physics.producers.mcgen.FluxFiles: ["gsimple_bnb_neutrino_icarus_dirt_*.root"]
9+
physics.producers.mcgen.FluxSearchPaths: "/cvmfs/sbn.osgstorage.org/pnfs/fnal.gov/usr/sbn/persistent/stash/physics-gputnam/icarus-bnb-dirt/"
10+
physics.producers.generator: @erase
11+
12+
# Reject neutrino interactions inside volDetEnclosure
13+
physics.filters.generator: @local::filtermctruthvolume
14+
15+
physics.simulate: [rns, mcgen, generator, beamgate]
16+
17+
outputs.rootoutput.outputCommands: [
18+
"keep *_*_*_*",
19+
"drop *_mcgen_*_*"
20+
]
21+
outputs.rootoutput.SelectEvents: ["simulate"]
Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
#include "simulation_genie_icarus_bnb_volDetEnclosure.fcl"
2+
3+
physics.producers.generator.TopVolume: "volWorld"
4+
physics.producers.generator.FiducialCut: "rockbox:(-378.49,-191.86,-904.950652270838)(+378.49,+144.96,+904.950652270838),1,800,0.00425,1.3,1"
5+
physics.producers.generator.FluxFiles: ["gsimple_bnb_neutrino_icarus_dirt_*.root"]
6+
physics.producers.generator.FluxSearchPaths: "/cvmfs/sbn.osgstorage.org/pnfs/fnal.gov/usr/sbn/persistent/stash/physics-gputnam/icarus-bnb-dirt/"
Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
#include "services_common_icarus.fcl"
2+
#include "icarus_siminfomixer.fcl"
3+
4+
services:
5+
{
6+
@table::icarus_common_services
7+
8+
scheduler: { defaultExceptions: false } # Make all uncaught exceptions fatal.
9+
# Load the service that manages root files for histograms.
10+
TFileService: { fileName: "siminfomixer.root" }
11+
12+
RandomNumberGenerator: {} #ART native random number generator
13+
}
14+
15+
process_name : SimInfoMixer #The process name must NOT contain any underscores
16+
17+
source:
18+
{
19+
module_type: RootInput
20+
saveMemoryObjectThreshold: 0
21+
maxEvents: -1
22+
}
23+
24+
outputs: {
25+
out: { module_type: RootOutput
26+
fileName: "%ifb_%tc_mixed.root"
27+
compressionLevel: 1
28+
dataTier: "generated"
29+
SelectEvents: ["mixer_path"]
30+
}
31+
}
32+
33+
physics: {
34+
35+
producers : {
36+
}
37+
38+
analyzers: {
39+
}
40+
41+
filters : {
42+
generator: @local::icarus_siminfomixer
43+
beamgate: @local::icarus_siminfomixer
44+
largeant: @local::icarus_siminfomixer
45+
ionization: @local::icarus_siminfomixer
46+
simplemerge: @local::icarus_siminfomixer
47+
genericcrt: @local::icarus_siminfomixer
48+
sedlite: @local::icarus_siminfomixer
49+
simdrift: @local::icarus_siminfomixer
50+
pdfastsim: @local::icarus_siminfomixer
51+
}
52+
53+
mixer_path : [ generator, beamgate, largeant, ionization, simplemerge, genericcrt, sedlite, simdrift, pdfastsim]
54+
trigger_paths : [ mixer_path ]
55+
56+
output : [ out ]
57+
end_paths: [ output ]
58+
59+
}
60+
61+
physics.filters.generator.MCTruthInputModuleLabels: ["generator::GenGenie"]
62+
physics.filters.generator.GTruthInputModuleLabels: ["generator::GenGenie"]
63+
physics.filters.generator.MCTruthGTruthAssnsInputModuleLabels: ["generator::GenGenie"]
64+
physics.filters.generator.MCFluxInputModuleLabels: ["generator::GenGenie"]
65+
physics.filters.generator.MCTruthMCFluxAssnsInputModuleLabels: ["generator::GenGenie"]
66+
physics.filters.generator.FillPOTInfo: true
67+
physics.filters.generator.POTSummaryTag: "potinevent:SubRunPOT:G4"
68+
69+
physics.filters.beamgate.BeamGateInputModuleLabels: ["beamgate::GenGenie"]
70+
71+
physics.filters.largeant.SimEnergyDepositInputModuleLabels: [
72+
"largeant:LArG4DetectorServicevolTPCPlaneV:G4",
73+
"largeant:LArG4DetectorServicevolTPC0:G4",
74+
"largeant:LArG4DetectorServicevolTPCPlaneU:G4",
75+
"largeant:LArG4DetectorServicevolTPCPlaneY:G4",
76+
"largeant:LArG4DetectorServicevolTPCActive:G4"
77+
]
78+
physics.filters.largeant.AuxDetHitInputModuleLabels: [
79+
"largeant:LArG4DetectorServicevolAuxDetSensitiveCERNbot:G4",
80+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut309:G4",
81+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOS:G4",
82+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut256:G4",
83+
"largeant:LArG4DetectorServicevolAuxDetSensitiveDC:G4",
84+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut508:G4",
85+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut497:G4",
86+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut325:G4",
87+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut485:G4",
88+
"largeant:LArG4DetectorServicevolAuxDetSensitiveCERNtop:G4",
89+
"largeant:LArG4DetectorServicevolAuxDetSensitiveMINOScut400:G4"
90+
]
91+
physics.filters.largeant.MCParticleInputModuleLabels: [
92+
"largeant::G4",
93+
"largeant:droppedMCParticles:G4"
94+
]
95+
physics.filters.largeant.MCTruthInputModuleLabels: ["generator::GenGenie"]
96+
physics.filters.largeant.MCTruthMCParticleAssnsInputModuleLabels: ["largeant::G4"]
97+
98+
physics.filters.ionization.SimEnergyDepositInputModuleLabels: [
99+
"ionization::G4",
100+
"ionization:priorSCE:G4"
101+
]
102+
103+
physics.filters.simplemerge.MCParticleInputModuleLabels: ["simplemerge::G4"]
104+
105+
physics.filters.genericcrt.AuxDetSimChannelInputModuleLabels: ["genericcrt::G4"]
106+
107+
physics.filters.sedlite.SimEnergyDepositLiteInputModuleLabels: ["sedlite::G4"]
108+
109+
physics.filters.simdrift.SimChannelInputModuleLabels: ["simdrift::G4"]
110+
111+
physics.filters.pdfastsim.SimPhotonsInputModuleLabels: ["pdfastsim::G4"]
Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
////////////////////////////////////////////////////////////////////////
2+
// FilterMCTruthVolume_module.cc
3+
//
4+
// art::EDFilter that selects events containing at least one MCTruth
5+
// neutrino interaction with a vertex __outside__ a user-specified box.
6+
// Produces filtered MCTruth, GTruth, and MCFlux collections with
7+
// their associations. Passes through POT summary at the SubRun level.
8+
////////////////////////////////////////////////////////////////////////
9+
10+
#include "art/Framework/Core/EDFilter.h"
11+
#include "art/Framework/Core/ModuleMacros.h"
12+
#include "art/Framework/Principal/Event.h"
13+
#include "art/Framework/Principal/Handle.h"
14+
#include "art/Framework/Principal/SubRun.h"
15+
#include "art/Persistency/Common/PtrMaker.h"
16+
#include "canvas/Persistency/Common/Assns.h"
17+
#include "canvas/Persistency/Common/FindOneP.h"
18+
#include "canvas/Utilities/InputTag.h"
19+
#include "fhiclcpp/ParameterSet.h"
20+
#include "messagefacility/MessageLogger/MessageLogger.h"
21+
22+
#include "larcoreobj/SummaryData/POTSummary.h"
23+
#include "nusimdata/SimulationBase/MCTruth.h"
24+
#include "nusimdata/SimulationBase/GTruth.h"
25+
#include "nusimdata/SimulationBase/MCFlux.h"
26+
27+
#include <array>
28+
#include <memory>
29+
#include <vector>
30+
31+
class FilterMCTruthVolume : public art::EDFilter {
32+
public:
33+
explicit FilterMCTruthVolume(fhicl::ParameterSet const& pset);
34+
35+
bool filter(art::Event& evt) override;
36+
bool endSubRun(art::SubRun& sr) override;
37+
38+
private:
39+
art::InputTag fMCTruthLabel;
40+
art::InputTag fPOTLabel;
41+
std::array<double, 3> fVolumeLow;
42+
std::array<double, 3> fVolumeHigh;
43+
44+
bool inVolume(double x, double y, double z) const;
45+
};
46+
47+
// ---------------------------------------------------------------------------
48+
FilterMCTruthVolume::FilterMCTruthVolume(fhicl::ParameterSet const& pset)
49+
: EDFilter{pset}
50+
, fMCTruthLabel{pset.get<art::InputTag>("MCTruthLabel", "generator")}
51+
, fPOTLabel{pset.get<art::InputTag>("POTLabel", "generator")}
52+
{
53+
auto lo = pset.get<std::vector<double>>("VolumeLow");
54+
auto hi = pset.get<std::vector<double>>("VolumeHigh");
55+
if (lo.size() != 3 || hi.size() != 3)
56+
throw art::Exception(art::errors::Configuration)
57+
<< "VolumeLow and VolumeHigh must each have exactly 3 elements (x, y, z).";
58+
std::copy(lo.begin(), lo.end(), fVolumeLow.begin());
59+
std::copy(hi.begin(), hi.end(), fVolumeHigh.begin());
60+
61+
produces<std::vector<simb::MCTruth>>();
62+
produces<std::vector<simb::GTruth>>();
63+
produces<std::vector<simb::MCFlux>>();
64+
produces<art::Assns<simb::MCTruth, simb::GTruth>>();
65+
produces<art::Assns<simb::MCTruth, simb::MCFlux>>();
66+
produces<sumdata::POTSummary, art::InSubRun>();
67+
}
68+
69+
// ---------------------------------------------------------------------------
70+
bool FilterMCTruthVolume::inVolume(double x, double y, double z) const
71+
{
72+
return x >= fVolumeLow[0] && x <= fVolumeHigh[0] &&
73+
y >= fVolumeLow[1] && y <= fVolumeHigh[1] &&
74+
z >= fVolumeLow[2] && z <= fVolumeHigh[2];
75+
}
76+
77+
// ---------------------------------------------------------------------------
78+
bool FilterMCTruthVolume::filter(art::Event& evt)
79+
{
80+
auto outMCTruth = std::make_unique<std::vector<simb::MCTruth>>();
81+
auto outGTruth = std::make_unique<std::vector<simb::GTruth>>();
82+
auto outMCFlux = std::make_unique<std::vector<simb::MCFlux>>();
83+
auto assnsTG = std::make_unique<art::Assns<simb::MCTruth, simb::GTruth>>();
84+
auto assnsTF = std::make_unique<art::Assns<simb::MCTruth, simb::MCFlux>>();
85+
86+
auto const& mctruthHandle = evt.getValidHandle<std::vector<simb::MCTruth>>(fMCTruthLabel);
87+
88+
art::FindOneP<simb::GTruth> findGTruth(mctruthHandle, evt, fMCTruthLabel);
89+
art::FindOneP<simb::MCFlux> findMCFlux(mctruthHandle, evt, fMCTruthLabel);
90+
91+
art::PtrMaker<simb::MCTruth> makeTruthPtr(evt);
92+
art::PtrMaker<simb::GTruth> makeGTruthPtr(evt);
93+
art::PtrMaker<simb::MCFlux> makeMCFluxPtr(evt);
94+
95+
for (size_t i = 0; i < mctruthHandle->size(); ++i) {
96+
auto const& mct = mctruthHandle->at(i);
97+
98+
if (!mct.NeutrinoSet()) continue;
99+
100+
double vx = mct.GetNeutrino().Nu().Vx();
101+
double vy = mct.GetNeutrino().Nu().Vy();
102+
double vz = mct.GetNeutrino().Nu().Vz();
103+
104+
if (inVolume(vx, vy, vz)) continue;
105+
106+
outMCTruth->push_back(mct);
107+
size_t idx = outMCTruth->size() - 1;
108+
109+
if (findGTruth.isValid()) {
110+
auto const& gtp = findGTruth.at(i);
111+
if (gtp.isNonnull()) {
112+
outGTruth->push_back(*gtp);
113+
assnsTG->addSingle(makeTruthPtr(idx), makeGTruthPtr(outGTruth->size() - 1));
114+
}
115+
}
116+
117+
if (findMCFlux.isValid()) {
118+
auto const& mfp = findMCFlux.at(i);
119+
if (mfp.isNonnull()) {
120+
outMCFlux->push_back(*mfp);
121+
assnsTF->addSingle(makeTruthPtr(idx), makeMCFluxPtr(outMCFlux->size() - 1));
122+
}
123+
}
124+
}
125+
126+
bool pass = !outMCTruth->empty();
127+
128+
evt.put(std::move(outMCTruth));
129+
evt.put(std::move(outGTruth));
130+
evt.put(std::move(outMCFlux));
131+
evt.put(std::move(assnsTG));
132+
evt.put(std::move(assnsTF));
133+
134+
return pass;
135+
}
136+
137+
// ---------------------------------------------------------------------------
138+
bool FilterMCTruthVolume::endSubRun(art::SubRun& sr)
139+
{
140+
auto const& potHandle = sr.getValidHandle<sumdata::POTSummary>(fPOTLabel);
141+
sr.put(std::make_unique<sumdata::POTSummary>(*potHandle), art::subRunFragment());
142+
return true;
143+
}
144+
145+
DEFINE_ART_MODULE(FilterMCTruthVolume)
Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
#include "art/Framework/Core/EDFilter.h"
2+
#include "art/Framework/Core/ModuleMacros.h"
3+
#include "art/Framework/Principal/Event.h"
4+
5+
#include "lardataobj/Simulation/SimEnergyDeposit.h"
6+
7+
namespace filt {
8+
9+
class SimEnergyDepFakeTriggerFilterICARUS : public art::EDFilter {
10+
public:
11+
explicit SimEnergyDepFakeTriggerFilterICARUS(fhicl::ParameterSet const& pset);
12+
virtual bool filter(art::Event& e) override;
13+
14+
private:
15+
const double fBeamTimeMin; // Minimum time of beam window [us]
16+
const double fBeamTimeMax; // Maximum time of beam window [us]
17+
const double fEnergyDeposit; // Minimum energy deposit in TPC for trigger [MeV]
18+
const double fMaxEnergyDeposit; // Maximum energy deposit in TPC for trigger [MeV]
19+
20+
const std::string fSimEnergyDepModuleName;
21+
};
22+
23+
SimEnergyDepFakeTriggerFilterICARUS::SimEnergyDepFakeTriggerFilterICARUS(fhicl::ParameterSet const& pset)
24+
: EDFilter(pset)
25+
, fBeamTimeMin(pset.get<double>("BeamTimeMin"))
26+
, fBeamTimeMax(pset.get<double>("BeamTimeMax"))
27+
, fEnergyDeposit(pset.get<double>("EnergyDeposit"))
28+
, fMaxEnergyDeposit(pset.get<double>("MaxEnergyDeposit", std::numeric_limits<double>::max()))
29+
, fSimEnergyDepModuleName(pset.get<std::string>("SimEnergyDepModuleName"))
30+
{
31+
}
32+
33+
bool SimEnergyDepFakeTriggerFilterICARUS::filter(art::Event& e)
34+
{
35+
const art::ValidHandle<std::vector<sim::SimEnergyDeposit>>&
36+
energyDeps(e.getValidHandle<std::vector<sim::SimEnergyDeposit>>(fSimEnergyDepModuleName));
37+
38+
double energy(0);
39+
40+
for (const sim::SimEnergyDeposit& energyDep : *energyDeps) {
41+
// Check particle time is within the beam time
42+
const double time(energyDep.Time() * 1e-3); // [ns] -> [us]
43+
if (time < fBeamTimeMin || time > fBeamTimeMax)
44+
continue;
45+
46+
// Add up the energy deposit inside the TPC
47+
energy += energyDep.Energy(); // [MeV]
48+
}
49+
50+
std::cout << "SAW E= " << energy << " inside window. Threshold is: " << fEnergyDeposit << std::endl;
51+
// If the energy deposit within the beam time is greater than some limit then trigger the event
52+
return (energy > fEnergyDeposit) & (energy < fMaxEnergyDeposit);
53+
}
54+
55+
DEFINE_ART_MODULE(SimEnergyDepFakeTriggerFilterICARUS)
56+
57+
}
58+
Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
BEGIN_PROLOG
2+
3+
filtermctruthvolume: {
4+
module_type: FilterMCTruthVolume
5+
MCTruthLabel: "mcgen"
6+
POTLabel: "mcgen"
7+
VolumeLow: [-571.55, -352.067, -1187.04]
8+
VolumeHigh: [+571.55, +619.033, +1566.2]
9+
}
10+
11+
END_PROLOG
Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
BEGIN_PROLOG
2+
3+
icarus_simenergydepfaketriggerfilter:
4+
{
5+
module_type: "SimEnergyDepFakeTriggerFilterICARUS"
6+
7+
BeamTimeMin: -0.2 # Minimum time of beam window [us]
8+
BeamTimeMax: 1.9 # Maximum time of beam window [us]
9+
EnergyDeposit: 5 # Minimum energy deposit in TPC for trigger [MeV]
10+
11+
# By default, take only the energy deposits within the TPC active volume
12+
SimEnergyDepModuleName: "largeant:LArG4DetectorServicevolTPCActive" # Name of SimEnergyDeposit producer module
13+
}
14+
15+
END_PROLOG

0 commit comments

Comments
 (0)