diff --git a/fcl/CMakeLists.txt b/fcl/CMakeLists.txt index f48ba48e4..d48921420 100644 --- a/fcl/CMakeLists.txt +++ b/fcl/CMakeLists.txt @@ -11,5 +11,7 @@ add_subdirectory(utilities) add_subdirectory(compatibility) add_subdirectory(numi-anaA) add_subdirectory(syst_variations) +add_subdirectory(decoder) +add_subdirectory(overlays) install_fhicl() diff --git a/fcl/caf/cafmakerjob_icarus_detsim2d_overlay.fcl b/fcl/caf/cafmakerjob_icarus_detsim2d_overlay.fcl new file mode 100644 index 000000000..171ce0b0b --- /dev/null +++ b/fcl/caf/cafmakerjob_icarus_detsim2d_overlay.fcl @@ -0,0 +1,7 @@ +#include "cafmakerjob_icarus.fcl" +#include "cafmaker_add_detsim2d_icarus.fcl" + +services.BackTrackerService.BackTracker.OverrideRealData: true +services.ParticleInventoryService.ParticleInventory.OverrideRealData: true + +physics.producers.cafmaker.OverrideRealData: true diff --git a/fcl/decoder/CMakeLists.txt b/fcl/decoder/CMakeLists.txt new file mode 100644 index 000000000..520ed808b --- /dev/null +++ b/fcl/decoder/CMakeLists.txt @@ -0,0 +1,8 @@ +# Install fcl files in /job subdirectory. + +install_fhicl() + +# Also put a copy in the source tree. + +FILE(GLOB fcl_files *.fcl) +install_source( EXTRAS ${fcl_files} ) diff --git a/fcl/decoder/decoder_data_icarus.fcl b/fcl/decoder/decoder_data_icarus.fcl new file mode 100644 index 000000000..852fc9855 --- /dev/null +++ b/fcl/decoder/decoder_data_icarus.fcl @@ -0,0 +1,56 @@ +### +## This fhicl file is used to run "stage0" processing specifically for the case where all +## TPC data is included in an artdaq data product from a single instance +## +#include "stage0_icarus_driver_common.fcl" + +process_name: decoder + +## Define the path we'll execute +physics.path: [ "filterdataintegrity", + "triggerconfig", + "daqTrigger", # Note sure if we actually need some of these trigger and PMT stuff other than daqPMT... Gianluca? + "pmtconfig", + "daqPMT", + "pmtconfigbaselines", + "pmtthr", # ??? + "pmtbaselines", + "daqCRT", + "crthit", + #"crtpmt", # Can we run this later: we probably want to, to have the PMT/CRT info for the nu + "daqTPCROI" + ] + +## boiler plate... +physics.outana: [ ] +physics.trigger_paths: [ path ] +physics.end_paths: [ outana, streamROOT ] + +# Drop the artdaq format files on output, +# Drop all output from the TPC decoder stage +# Drop all output from the 1D deconvolution stage +# Drop the recob::Wire output from the roifinder (but keep the ChannelROIs) +# Drop the reconstructed optical hits before the timing correction +# Drop the PMT waveforms (will keep the ones from daqPMTonbeam) +# Keep the Trigger fragment +outputs.rootOutput.outputCommands: [ + "keep *_*_*_*", + "drop *_*_*_DAQ*", + "drop *_ophituncorrected_*_*", + #"drop raw::OpDetWaveforms_daqPMT__*", + "drop *_decon1droi_*_*", + "drop recob::Wire*_roifinder_*_*", + "keep *_daq_ICARUSTriggerV*_*", + "drop recob::ChannelROIs_daqTPCROI_*_*", #ChannelROIs_daqTPCROI_PHYSCRATEDATATPCWE_decoder decoder..... | daqTPCROI......... | PHYSCRATEDATATPCWW + "drop raw::RawDigits_daqTPCROI_PHYSCRATEDATATPCWE_decoder", #RawDigits_daqTPCROI_PHYSCRATEDATATPCWW_decoder.obj + "drop raw::RawDigits_daqTPCROI_PHYSCRATEDATATPCWW_decoder", + "drop raw::RawDigits_daqTPCROI_PHYSCRATEDATATPCEW_decoder", + "drop raw::RawDigits_daqTPCROI_PHYSCRATEDATATPCEE_decoder", + "drop artdaq::Fragments_daq_*_ICARUSReprocessRaw" + ] + +## Modify the event selection for the purity analyzers +physics.analyzers.purityinfoana0.SelectEvents: [ path ] +physics.analyzers.purityinfoana1.SelectEvents: [ path ] + +physics.producers.daqTPCROI.OutputRawWaveform: true diff --git a/fcl/decoder/decoder_run2_icarus.fcl b/fcl/decoder/decoder_run2_icarus.fcl new file mode 100644 index 000000000..4a92eea4f --- /dev/null +++ b/fcl/decoder/decoder_run2_icarus.fcl @@ -0,0 +1,48 @@ +### +## This fhicl file is used to run "stage0" processing specifically for the case where all +## TPC data is included in an artdaq data product from a single instance +## +#include "stage0_icarus_driver_common.fcl" + +process_name: stage0 + +## Define the path we'll execute +physics.path: [ "filterdataintegrity", + "triggerconfig", + "daqTrigger", # Note sure if we actually need some of these trigger and PMT stuff other than daqPMT... Gianluca? + "pmtconfig", + "daqPMT", + "pmtconfigbaselines", + "pmtthr", # ??? + "pmtbaselines", + "daqCRT", + "crthit", + #"crtpmt", # Can we run this later: we probably want to, to have the PMT/CRT info for the nu + "daqTPCROI" ] + +## boiler plate... +physics.outana: [ ] +physics.trigger_paths: [ path ] +physics.end_paths: [ outana, streamROOT ] + +# Drop the artdaq format files on output, +# Drop all output from the TPC decoder stage +# Drop all output from the 1D deconvolution stage +# Drop the recob::Wire output from the roifinder (but keep the ChannelROIs) +# Drop the reconstructed optical hits before the timing correction +# Drop the PMT waveforms (will keep the ones from daqPMTonbeam) +# Keep the Trigger fragment +outputs.rootOutput.outputCommands: [ + "keep *_*_*_*", + "drop *_*_*_DAQ*", + "drop *_ophituncorrected_*_*", + #"drop raw::OpDetWaveforms_daqPMT__*", + "drop *_decon1droi_*_*", + "drop recob::Wire*_roifinder_*_*", + "keep *_daq_ICARUSTriggerV*_*"] + +## Modify the event selection for the purity analyzers +physics.analyzers.purityinfoana0.SelectEvents: [ path ] +physics.analyzers.purityinfoana1.SelectEvents: [ path ] + +physics.producers.daqTPCROI.OutputRawWaveform: true diff --git a/fcl/detsim/detsim_2d_icarus_fitFR_nonoise.fcl b/fcl/detsim/detsim_2d_icarus_fitFR_nonoise.fcl new file mode 100644 index 000000000..efc73667f --- /dev/null +++ b/fcl/detsim/detsim_2d_icarus_fitFR_nonoise.fcl @@ -0,0 +1,5 @@ +#include "detsim_2d_icarus_fitFR.fcl" + + +physics.producers.daq.wcls_main.structs.coh_noise_scale: 0.0 +physics.producers.daq.wcls_main.structs.int_noise_scale: 0.0 diff --git a/fcl/detsim/detsim_2d_icarus_fitFR_overlay.fcl b/fcl/detsim/detsim_2d_icarus_fitFR_overlay.fcl new file mode 100644 index 000000000..b66d01caa --- /dev/null +++ b/fcl/detsim/detsim_2d_icarus_fitFR_overlay.fcl @@ -0,0 +1,9 @@ +#include "detsim_2d_icarus_fitFR.fcl" + + +physics.producers.daq.wcls_main.structs.coh_noise_scale: 0.0 +physics.producers.daq.wcls_main.structs.int_noise_scale: 0.0 +physics.producers.daq.wcls_main.structs.include_noise: false +physics.producers.daq.wcls_main.structs.overlay_drifter: true +physics.producers.daq.wcls_main.plugins: [@sequence::physics.producers.daq.wcls_main.plugins, "WireCellICARUSDrifter", "WireCellHio"] +physics.producers.daq.wcls_main.inputers: ["wclsICARUSDrifter", "wclsSimDepoSetSource:electron"] diff --git a/fcl/overlays/CMakeLists.txt b/fcl/overlays/CMakeLists.txt new file mode 100644 index 000000000..520ed808b --- /dev/null +++ b/fcl/overlays/CMakeLists.txt @@ -0,0 +1,8 @@ +# Install fcl files in /job subdirectory. + +install_fhicl() + +# Also put a copy in the source tree. + +FILE(GLOB fcl_files *.fcl) +install_source( EXTRAS ${fcl_files} ) diff --git a/fcl/overlays/icarus_siminfomixer.fcl b/fcl/overlays/icarus_siminfomixer.fcl new file mode 100644 index 000000000..8f381e326 --- /dev/null +++ b/fcl/overlays/icarus_siminfomixer.fcl @@ -0,0 +1,33 @@ +BEGIN_PROLOG + +icarus_siminfomixer: +{ + module_type : SimInfoOverlayFilter + SimInputFileNames : [ "genfile.root.local" ] + AuxDetSimChannelInputModuleLabels: [] + SimEnergyDepositInputModuleLabels: [] + MCParticleInputModuleLabels: [] + MCTruthGTruthAssnsInputModuleLabels: [] + MCTruthInputModuleLabels: [] + SimChannelInputModuleLabels: [] + SimPhotonsInputModuleLabels: [] + MCTruthMCParticleAssnsMap: [] + MCTruthMCParticleAssnsInputModuleLabels: [] + BeamGateInputModuleLabels: [] + MCParticleLiteInputModuleLabels: [] + MCFluxInputModuleLabels: [] + MCTruthMCFluxAssnsInputModuleLabels: [] + + Verbosity: -1 + FillPOTInfo: false + POTSummaryTag: "generator::GenieGenFiltered" + detail : {} +} + +icarus_subrunpotinevent: +{ + module_type : SubRunPOTInEvent + InputTag: "generator:SubRunPOT" +} + +END_PROLOG \ No newline at end of file diff --git a/fcl/overlays/intime_gen_overlay_SimInfoMixer1.fcl b/fcl/overlays/intime_gen_overlay_SimInfoMixer1.fcl new file mode 100644 index 000000000..25a7ee569 --- /dev/null +++ b/fcl/overlays/intime_gen_overlay_SimInfoMixer1.fcl @@ -0,0 +1,75 @@ +#include "icarus_siminfomixer.fcl" +#include "local_gen_include.fcl" + +//======================================== + +services: +{ + scheduler: { defaultExceptions: false } # Make all uncaught exceptions fatal. + # Load the service that manages root files for histograms. + TFileService: { fileName: "gen_hist.root" } + + RandomNumberGenerator: {} #ART native random number generator + FileCatalogMetadata: @local::art_file_catalog_mc + #@table::microboone_g4_services +} + +process_name : SimInfoMixer1 #The process name must NOT contain any underscores + +source: +{ + module_type: RootInput + saveMemoryObjectThreshold: 0 + maxEvents: -1 +} + +outputs: { + out: { module_type: RootOutput + fileName: "%ifb_%tc_simmxd1.root" + compressionLevel: 1 + dataTier: "generated" + SelectEvents: ["mixer_path"] + } +} + +physics: { + + producers : { + } + + analyzers: { + } + + filters : { + generator: @local::icarus_siminfomixer + beamgate: @local::icarus_siminfomixer + GenInTimeSorter: @local::icarus_siminfomixer + larg4intime: @local::icarus_siminfomixer + } + + mixer_path : [ generator, beamgate, GenInTimeSorter, larg4intime ] + trigger_paths : [ mixer_path ] + + output : [ out ] + end_paths: [ output ] + +} + +##:: 2:3:1 +##PROCESS NAME...... | MODULE LABEL... | PRODUCT INSTANCE NAME.... | DATA PRODUCT TYPE.................................................... | .SIZE +physics.filters.generator.SimInputFileNames : [ "genfile.root.local" ] + +physics.filters.larg4intime.AuxDetSimChannelInputModuleLabels: ["larg4intime::CosmicsCorsikaPGenAndG4InTime"] +physics.filters.larg4intime.MCParticleLiteInputModuleLabels: ["larg4intime::CosmicsCorsikaPGenAndG4InTime"] + +physics.filters.larg4intime.MCParticleInputModuleLabels: ["larg4intime::CosmicsCorsikaPGenAndG4InTime"] +physics.filters.larg4intime.MCTruthMCParticleAssnsInputModuleLabels: ["larg4intime::CosmicsCorsikaPGenAndG4InTime"] +physics.filters.larg4intime.SimEnergyDepositInputModuleLabels: ["larg4intime:TPCActive:CosmicsCorsikaPGenAndG4InTime"] + +physics.filters.beamgate.BeamGateInputModuleLabels: ["beamgate::CosmicsCorsikaPGenAndG4InTime"] +physics.filters.GenInTimeSorter.MCTruthInputModuleLabels: ["GenInTimeSorter:intime:CosmicsCorsikaPGenAndG4InTime"] #std::vector intime +physics.filters.larg4intime.SimChannelInputModuleLabels: ["larg4intime::CosmicsCorsikaPGenAndG4InTime"] +physics.filters.generator.MCTruthInputModuleLabels: ["GenInTimeSorter:intime:CosmicsCorsikaPGenAndG4InTime"] #std::vector #no label +physics.filters.generator.MCTruthInputModuleLabels: ["GenInTimeSorter::CosmicsCorsikaPGenAndG4InTime"] #std::vector #no label +physics.filters.generator.MCTruthInputModuleLabels: ["generator::CosmicsCorsikaPGenAndG4InTime"] +physics.filters.larg4intime.SimPhotonsInputModuleLabels: ["larg4intime::CosmicsCorsikaPGenAndG4InTime"] diff --git a/fcl/overlays/intime_gen_overlay_SimInfoMixer2.fcl b/fcl/overlays/intime_gen_overlay_SimInfoMixer2.fcl new file mode 100644 index 000000000..bf5b1b743 --- /dev/null +++ b/fcl/overlays/intime_gen_overlay_SimInfoMixer2.fcl @@ -0,0 +1,58 @@ +#include "icarus_siminfomixer.fcl" +#include "local_gen_include.fcl" + +//======================================== + +services: +{ + scheduler: { defaultExceptions: false } # Make all uncaught exceptions fatal. + # Load the service that manages root files for histograms. + TFileService: { fileName: "gen_hist.root" } + RandomNumberGenerator: {} #ART native random number generator + FileCatalogMetadata: @local::art_file_catalog_mc +} + +process_name : SimInfoMixer2 #The process name must NOT contain any underscores + +source: +{ + module_type: RootInput + saveMemoryObjectThreshold: 0 + maxEvents: -1 +} + +outputs: { + out: { module_type: RootOutput + fileName: "%ifb_%tc_simmxd2.root" + compressionLevel: 1 + dataTier: "generated" + SelectEvents: ["mixer_path"] + } +} + +physics: { + + producers : { + } + + analyzers: { + } + + filters : { + GenInTimeSorter: @local::icarus_siminfomixer + larg4intime: @local::icarus_siminfomixer + } + + mixer_path : [ GenInTimeSorter, larg4intime ] + trigger_paths : [ mixer_path ] + + output : [ out ] + end_paths: [ output ] + +} + +##:: +physics.filters.larg4intime.SimInputFileNames : [ "genfile.root.local" ] + +physics.filters.larg4intime.SimEnergyDepositInputModuleLabels: ["larg4intime:Other:CosmicsCorsikaPGenAndG4InTime"] #Other +physics.filters.GenInTimeSorter.MCTruthInputModuleLabels: ["GenInTimeSorter:outtime:CosmicsCorsikaPGenAndG4InTime"] #std::vector outtime diff --git a/fcl/overlays/overlay_waveforms.fcl b/fcl/overlays/overlay_waveforms.fcl new file mode 100644 index 000000000..52d31b94c --- /dev/null +++ b/fcl/overlays/overlay_waveforms.fcl @@ -0,0 +1,105 @@ +### +## This fhicl file is used to run "stage0" processing specifically for the case where all +## TPC data is included in an artdaq data product from a single instance +## +#include "stage0_icarus_mc_defs.fcl" +#include "stage0_icarus_driver_common.fcl" + +process_name: Overlay + +## Add the MC module to the list of producers +physics.producers: { + crthit: @local::standard_crtsimhitproducer + overlayTPCRawEE: + { + module_type: "OverlayProducts" + TPCOverlayRaw: true + TPCOverlayHits: false + TPCHitsWireAssn: false + TPCOverlayROI: false + PMTOverlayRaw: false + PMTOverlayHits: false + CRTOverlayHits: false + TPCRawInputLabels: ["daqTPCROI:PHYSCRATEDATATPCEERAW", "daq:TPCEE"] + } + overlayTPCRawEW: + { + module_type: "OverlayProducts" + TPCOverlayRaw: true + TPCOverlayHits: false + TPCHitsWireAssn: false + TPCOverlayROI: false + PMTOverlayRaw: false + PMTOverlayHits: false + CRTOverlayHits: false + TPCRawInputLabels: ["daqTPCROI:PHYSCRATEDATATPCEWRAW", "daq:TPCEW"] + } + overlayTPCRawWE: + { + module_type: "OverlayProducts" + TPCOverlayRaw: true + TPCOverlayHits: false + TPCHitsWireAssn: false + TPCOverlayROI: false + PMTOverlayRaw: false + PMTOverlayHits: false + CRTOverlayHits: false + TPCRawInputLabels: ["daqTPCROI:PHYSCRATEDATATPCWERAW", "daq:TPCWE"] + } + overlayTPCRawWW: + { + module_type: "OverlayProducts" + TPCOverlayRaw: true + TPCOverlayHits: false + TPCHitsWireAssn: false + TPCOverlayROI: false + PMTOverlayRaw: false + PMTOverlayHits: false + CRTOverlayHits: false + TPCRawInputLabels: ["daqTPCROI:PHYSCRATEDATATPCWWRAW", "daq:TPCWW"] + } + + overlayOpWaveforms: + { + module_type: "OverlayProducts" + TPCOverlayRaw: false + TPCOverlayHits: false + TPCHitsWireAssn: false + TPCOverlayROI: false + PMTOverlayRaw: true + PMTOverlayHits: false + CRTOverlayHits: false + PMTWaveDataLabel: "daqPMT" + PMTWaveSimLabel: "opdaq" + PMTWaveBaseLabel: "pmtbaselines" + } + + overlayCRTHit: + { + module_type: "OverlayProducts" + TPCOverlayRaw: false + TPCOverlayHits: false + TPCHitsWireAssn: false + TPCOverlayROI: false + PMTOverlayRaw: false + PMTOverlayHits: false + CRTOverlayHits: true + CRTHitInputLabels: ["crthit::decoder","crthit::Overlay"] + } + } + +## Use the following to run the full defined stage0 set of modules +physics.reco: [ crthit, overlayTPCRawWW, overlayTPCRawWE, overlayTPCRawEW, overlayTPCRawEE, overlayOpWaveforms, overlayCRTHit] + +## boiler plate... +physics.outana: [ ] +physics.trigger_paths: [ reco ] +physics.end_paths: [ outana, streamROOT ] + +## Note that for output we hijack the "rootOutput" definition (but change the naming convention to make more generic for MC) +outputs.rootOutput.fileName: "%ifb_%tc-%p.root" +outputs.rootOutput.dataTier: "reconstructed" +outputs.rootOutput.SelectEvents: ["reco"] + +# Drop the artdaq format files on output +outputs.rootOutput.outputCommands: ["keep *_*_*_*","drop *_gaushitTPC*_*_*","drop *_ophit*_*_*","drop *_opflash*_*_*","drop *_crthit*_*_*","drop *_crttrack*_*_*"] #TODO: drop raw waveforms that arent from this module. Keep wf from this module to feed to later stages diff --git a/fcl/overlays/simulation_genie_icarus_bnb_overlays_volDetEnclosure.fcl b/fcl/overlays/simulation_genie_icarus_bnb_overlays_volDetEnclosure.fcl new file mode 100644 index 000000000..faed055cd --- /dev/null +++ b/fcl/overlays/simulation_genie_icarus_bnb_overlays_volDetEnclosure.fcl @@ -0,0 +1,24 @@ +# File: simulation_genie_icarus_bnb_volDetEnclosure.fcl +# Purpose: generate neutrino interactions from BNB in the whole DetEnclosure volume +# Date: January 11, 2021 +# Author: Gianluca Petrillo (petrillo@slac.stanford.edu) +# +# This configuration is based on `simulation_genie_icarus_bnb.fcl`, +# the only difference being that the generation volume is explicitly +# specified to be the detector enclosure, that is a box including +# all cryostats. +# + +#include "simulation_genie_icarus_bnb.fcl" + +physics.producers.generator.TopVolume: "volDetEnclosure" + +outputs.rootoutput.fileName: "simulation_genie_icarus_bnb_volDetEnclosure_%tc-%p.root" +services.TFileService.fileName: "Supplemental-simulation_genie_icarus_bnb_volDetEnclosure_%tc-%p.root" + +source: { + module_type: RootInput + maxEvents: -1 # Number of events to create + saveMemoryObjectThreshold: 0 +} # source + diff --git a/fcl/overlays/stage0_run2_icarus_overlay.fcl b/fcl/overlays/stage0_run2_icarus_overlay.fcl new file mode 100644 index 000000000..b9538a555 --- /dev/null +++ b/fcl/overlays/stage0_run2_icarus_overlay.fcl @@ -0,0 +1,64 @@ +### +## This fhicl file is used to run "stage0" processing specifically for the case where all +## TPC data is included in an artdaq data product from a single instance +## +#include "stage0_icarus_mc_defs.fcl" +#include "stage0_icarus_driver_common.fcl" + +process_name: MCstage0 + +## Add the MC module to the list of producers +physics.producers: { @table::icarus_stage0_producers + @table::icarus_stage0_mc_producers + } + +## Use the following to run the full defined stage0 set of modules +physics.path: [ @sequence::icarus_stage0_mc_PMT, + MCDecodeTPCROI, + @sequence::icarus_stage0_multiTPC, + @sequence::icarus_stage0_mc_crt + ] + +## boiler plate... +physics.outana: [ purityinfoana0, purityinfoana1 ] +physics.trigger_paths: [ path ] +physics.end_paths: [ outana, streamROOT ] + +# Drop data products that are no longer needed, but make sure to keep important items! +# For example, we need to drop the RawDigits from the detector simulation stage but want to keep the SimChannel info from WireCell... +#outputs.rootOutput.outputCommands: ["keep *_*_*_*", "drop *_daq*_*_*", "drop *_MCDecodeTPCROI_*_*", "drop *_decon1droi_*_*", "drop recob::Wire*_roifinder_*_*", "keep *_daq_simpleSC_*"] +outputs.rootOutput.outputCommands: ["keep *_*_*_*", + "drop *_daq*_*_*", + "drop *_MCDecodeTPCROI_*_*", + "drop *_decon1droi_*_*", + "drop recob::Wire*_roifinder_*_*", + "keep *_daq_simpleSC_*", + "drop raw::RawDigits_overlayTPCRaw*_*_Overlay", + "drop sim::SimEnergyDeposits_largeant_TPCActive_G4", + "drop sim::SimEnergyDeposits_shifted_*_DetSim", + "drop raw::OpDetWaveforms_overlayOpWaveforms_*_Overlay" + ] +# Set the expected input for ophit +physics.producers.ophit.InputModule: "overlayOpWaveforms" + +# Note the default assumption is that our detector simulation input will come from the WireCell 2D drift simulation, a la 'daq' +# If we are running the 1D drift simulation we need to switch to using: +# `physics.producers.MCDecodeTPCROI.FragmentsLabelVec: ["daq3:PHYSCRATEDATATPCWW","daq2:PHYSCRATEDATATPCWE","daq1:PHYSCRATEDATATPCEW","daq0:PHYSCRATEDATATPCEE"]` +# + +physics.producers.MCDecodeTPCROI.FragmentsLabelVec: ["overlayTPCRawWW", "overlayTPCRawWE", "overlayTPCRawEW", "overlayTPCRawEE"] +physics.producers.MCDecodeTPCROI.OutInstanceLabelVec: ["PHYSCRATEDATATPCWW", "PHYSCRATEDATATPCWE", "PHYSCRATEDATATPCEW", "PHYSCRATEDATATPCEE"] +physics.producers.decon1droi.RawDigitLabelVec: ["MCDecodeTPCROI:PHYSCRATEDATATPCWW","MCDecodeTPCROI:PHYSCRATEDATATPCWE","MCDecodeTPCROI:PHYSCRATEDATATPCEW","MCDecodeTPCROI:PHYSCRATEDATATPCEE"] +physics.producers.roifinder.WireModuleLabelVec: ["decon1droi:PHYSCRATEDATATPCWW","decon1droi:PHYSCRATEDATATPCWE","decon1droi:PHYSCRATEDATATPCEW","decon1droi:PHYSCRATEDATATPCEE"] + +physics.producers.decon2droiEE.wcls_main.params.raw_input_label: "MCDecodeTPCROI:PHYSCRATEDATATPCEE" +physics.producers.decon2droiEW.wcls_main.params.raw_input_label: "MCDecodeTPCROI:PHYSCRATEDATATPCEW" +physics.producers.decon2droiWE.wcls_main.params.raw_input_label: "MCDecodeTPCROI:PHYSCRATEDATATPCWE" +physics.producers.decon2droiWW.wcls_main.params.raw_input_label: "MCDecodeTPCROI:PHYSCRATEDATATPCWW" + +## Need overrides for the purity monitor +physics.analyzers.purityinfoana0.SelectEvents: [ path ] +physics.analyzers.purityinfoana1.SelectEvents: [ path ] +physics.producers.purityana0.RawModuleLabel: ["MCDecodeTPCROI:PHYSCRATEDATATPCWW","MCDecodeTPCROI:PHYSCRATEDATATPCWE","MCDecodeTPCROI:PHYSCRATEDATATPCEW","MCDecodeTPCROI:PHYSCRATEDATATPCEE"] +physics.producers.purityana1.RawModuleLabel: ["MCDecodeTPCROI:PHYSCRATEDATATPCWW","MCDecodeTPCROI:PHYSCRATEDATATPCWE","MCDecodeTPCROI:PHYSCRATEDATATPCEW","MCDecodeTPCROI:PHYSCRATEDATATPCEE"] + diff --git a/fcl/overlays/stage0_run2_icarus_overlay_wf.fcl b/fcl/overlays/stage0_run2_icarus_overlay_wf.fcl new file mode 100644 index 000000000..15de362a0 --- /dev/null +++ b/fcl/overlays/stage0_run2_icarus_overlay_wf.fcl @@ -0,0 +1,10 @@ +#include "stage0_run2_icarus_overlay.fcl" + +outputs.rootOutput.outputCommands: [ + "keep *_*_*_*", + "drop *_daq*_*_*", + "drop *_MCDecodeTPCROI_*_*", + "drop *_decon1droi_*_*", + "drop recob::Wire*_roifinder_*_*", + "keep *_daq_simpleSC_*" + ] \ No newline at end of file diff --git a/fcl/overlays/stage1_run2_1d_icarus_overlay.fcl b/fcl/overlays/stage1_run2_1d_icarus_overlay.fcl new file mode 100644 index 000000000..ae4ee90f0 --- /dev/null +++ b/fcl/overlays/stage1_run2_1d_icarus_overlay.fcl @@ -0,0 +1,4 @@ +#include "stage1_run2_icarus_overlay.fcl" + +physics.producers.cluster3DCryoW.Hit3DBuilderAlg.HitFinderTagVec: ["gaushit1dTPCWW", "gaushit1dTPCWE"] +physics.producers.cluster3DCryoE.Hit3DBuilderAlg.HitFinderTagVec: ["gaushit1dTPCEW", "gaushit1dTPCEE"] \ No newline at end of file diff --git a/fcl/overlays/stage1_run2_icarus_overlay.fcl b/fcl/overlays/stage1_run2_icarus_overlay.fcl new file mode 100644 index 000000000..1e24999ba --- /dev/null +++ b/fcl/overlays/stage1_run2_icarus_overlay.fcl @@ -0,0 +1,101 @@ + +#include "mchitmodules.fcl" +#include "mctrutht0matching.fcl" +#include "backtrackerservice.fcl" +#include "particleinventoryservice.fcl" +#include "stage1_icarus_driver_common.fcl" + +process_name: MCstage1 + +## Add the MC module to the list of producers +physics.producers: { + @table::icarus_stage1_producers + + mcophit: @local::ICARUSMCOpHit + mcopflashTPC0: @local::ICARUSMCOpFlashTPC0 + mcopflashTPC1: @local::ICARUSMCOpFlashTPC1 + mcopflashTPC2: @local::ICARUSMCOpFlashTPC2 + mcopflashTPC3: @local::ICARUSMCOpFlashTPC3 + + cheatopflashTPC0: @local::ICARUSCheatOpFlashTPC0 + cheatopflashTPC1: @local::ICARUSCheatOpFlashTPC1 + cheatopflashTPC2: @local::ICARUSCheatOpFlashTPC2 + cheatopflashTPC3: @local::ICARUSCheatOpFlashTPC3 + + ### mc producers + mchitfinder: @local::standard_mchitfinder + #mcassociationsGausCryoE: @local::standard_mcparticlehitmatching + #mcassociationsGausCryoW: @local::standard_mcparticlehitmatching +} + +physics.reco: [ + @sequence::icarus_reco_Gauss_CryoE , + @sequence::icarus_reco_Gauss_CryoW , + @sequence::icarus_reco_fm, + caloskimCalorimetryCryoE, caloskimCalorimetryCryoW + #mcassociationsGausCryoE, mcassociationsGausCryoW + ] + +# Turn on truth-info for track skimmer +physics.analyzers.caloskimE.G4producer: "largeant" +physics.analyzers.caloskimE.SimChannelproducer: "daq:simpleSC" +physics.analyzers.caloskimE.RawDigitproducers: ["MCDecodeTPCROI:PHYSCRATEDATATPCEW", "MCDecodeTPCROI:PHYSCRATEDATATPCEE"] +physics.analyzers.caloskimW.G4producer: "largeant" +physics.analyzers.caloskimW.SimChannelproducer: "daq:simpleSC" +physics.analyzers.caloskimW.RawDigitproducers: ["MCDecodeTPCROI:PHYSCRATEDATATPCWW", "MCDecodeTPCROI:PHYSCRATEDATATPCWE"] + +physics.outana: [ @sequence::icarus_analysis_modules ] +physics.trigger_paths: [ reco ] +physics.end_paths: [ outana, stream1 ] +outputs.out1.fileName: "%ifb_%tc-%p.root" +outputs.out1.dataTier: "reconstructed" +outputs.out1.outputCommands: [ + "keep *_*_*_*", + "drop *_caloskimCalorimetryCryoE_*_*", + "drop *_caloskimCalorimetryCryoW_*_*" +] + +# Disabled Space-Charge service for calorimetry +services.SpaceChargeService: { + EnableCalEfieldSCE: false + EnableCalSpatialSCE: false + EnableCorrSCE: false + EnableSimEfieldSCE: false + EnableSimSpatialSCE: false + InputFilename: "SCEoffsets/SCEoffsets_ICARUS_E500_voxelTH3.root" + RepresentationType: "Voxelized_TH3" + service_provider: "SpaceChargeServiceICARUS" +} + +services.BackTrackerService: @local::standard_backtrackerservice +# In the 2D-detsim, SimChannel objects are made by the WireCell +# drift simulation (daq), not LArG4 (largeant). Thus, we need +# to overwrite the truth matching labels in the calibration ntuple maker +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "daq:simpleSC" +services.ParticleInventoryService: @local::standard_particleinventoryservice + +services.BackTrackerService.BackTracker.OverrideRealData: true +services.ParticleInventoryService.ParticleInventory.OverrideRealData: true + + + +#physics.producers.mcassociationsGausCryoE.HitParticleAssociations.HitModuleLabelVec: ["cluster3DCryoE"] +#physics.producers.mcassociationsGausCryoW.HitParticleAssociations.HitModuleLabelVec: ["cluster3DCryoW"] + +services.message.destinations : +{ + STDCOUT: + { + type: "cout" #tells the message service to output this destination to cout + threshold: "WARNING" #tells the message service that this destination applies to WARNING and higher level messages + categories: + { + Cluster3DICARUS: + { + limit: -1 + reportEvery: 1 + } + } + } +} + diff --git a/icaruscode/CMakeLists.txt b/icaruscode/CMakeLists.txt index a87cbffd7..606f87f42 100644 --- a/icaruscode/CMakeLists.txt +++ b/icaruscode/CMakeLists.txt @@ -14,4 +14,5 @@ add_subdirectory(PMT) add_subdirectory(RecoUtils) add_subdirectory(TPC) add_subdirectory(Utilities) +add_subdirectory(Overlays) #add_subdirectory(Supera) diff --git a/icaruscode/Decode/CMakeLists.txt b/icaruscode/Decode/CMakeLists.txt index c64463b60..906e433ba 100644 --- a/icaruscode/Decode/CMakeLists.txt +++ b/icaruscode/Decode/CMakeLists.txt @@ -116,6 +116,22 @@ cet_build_plugin(CopyDaqICARUSTPC art::module LIBRARIES cetlib_except::cetlib_except ) +cet_build_plugin(OverlayProducts art::module LIBRARIES + icaruscode::Timing + sbnobj::Common_PMT_Data + larcore::Geometry_Geometry_service + art_root_io::TFileService_service + art_root_io::tfile_support + art::Framework_Services_Registry + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib_except::cetlib_except + ROOT::Tree + lardata::Utilities + lardataobj::RecoBase + lardata::ArtDataHelper + ) + install_headers() install_fhicl() install_source() diff --git a/icaruscode/Decode/OverlayProducts_module.cc b/icaruscode/Decode/OverlayProducts_module.cc new file mode 100644 index 000000000..c614a2eba --- /dev/null +++ b/icaruscode/Decode/OverlayProducts_module.cc @@ -0,0 +1,740 @@ +//////////////////////////////////////////////////////////////////////// +// Class: OverlayProducts +// Plugin Type: producer (Unknown Unknown) +// File: OverlayProducts_module.cc +// +// Generated at Fri Sep 23 17:04:48 2022 by Bruce Howard using cetskelgen +// from version . +//////////////////////////////////////////////////////////////////////// + +// NB: See https://github.com/SBNSoftware/icaruscode/blob/4ea3861350628402b7edccfe46efdd403e4c0432/icaruscode/TPC/SignalProcessing/RecoWire/ROIConverter_module.cc + +#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 "canvas/Persistency/Common/FindMany.h" +#include "canvas/Persistency/Common/FindManyP.h" + +#include "cetlib_except/exception.h" + +#include "lardataobj/RawData/RawDigit.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardata/ArtDataHelper/HitCreator.h" +//#include "icaruscode/IcarusObj/ChannelROI.h" +#include "sbnobj/ICARUS/TPC/ChannelROI.h" +#include "icaruscode/IcarusObj/classes.h" + +#include "lardataobj/RawData/OpDetWaveform.h" +#include "lardataobj/RecoBase/OpHit.h" +#include "sbnobj/ICARUS/PMT/Data/WaveformBaseline.h" + +// PROBABLY TEMPORARY +//#include "art/Framework/Services/Registry/ServiceHandle.h" +//#include "larcore/CoreUtils/ServiceUtil.h" +//#include "larcore/Geometry/Geometry.h" +///////////////////// + +#include "sbnobj/Common/CRT/CRTHit.hh" + +#include "art/Persistency/Common/PtrMaker.h" + +#include +#include +#include + +class OverlayProducts; + + +class OverlayProducts : public art::EDProducer { +public: + explicit OverlayProducts(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. + OverlayProducts(OverlayProducts const&) = delete; + OverlayProducts(OverlayProducts&&) = delete; + OverlayProducts& operator=(OverlayProducts const&) = delete; + OverlayProducts& operator=(OverlayProducts&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + + struct opWaveformOverlay { + raw::TimeStamp_t timestamp; + std::vector< raw::ADC_Count_t > wvfm; + float baseline; // TODO: does SBND have/plan to have such a product? + double deltaT = 0.002; // TODO: automate? SBND? + + raw::TimeStamp_t endTime() const { return timestamp+(deltaT*wvfm.size()); } + std::pair< raw::TimeStamp_t, raw::TimeStamp_t > timeBin( unsigned int binNumber ) { return std::make_pair(timestamp+(deltaT*binNumber), + timestamp+(deltaT*(binNumber+1))); } + }; + +private: + + // Declare member data here. + // TPC + bool fTPCOverlayRaw; //< Set true if you should overlay at the RawDigit stage + bool fTPCOverlayROI; //< Set true if you should overlay at the recob::ChannelROI stage + bool fTPCOverlayHits; //< Set true if you should overlay at the recob::Hit stage + bool fTPCHitsWireAssn; //< Set true if you want to make the hit<->wire association in HitCollectionCreator + std::vector< art::InputTag > fTPCRawInputLabels; + art::InputTag fTPCRawOutputLabel; + std::vector< art::InputTag > fTPCROIInputLabels; + std::vector< art::InputTag > fTPCHitInputLabels; + std::string fTPCHitCreatorInstanceName; + + // PMT + bool fPMTOverlayRaw; + bool fPMTOverlayHits; + art::InputTag fPMTWaveDataLabel; + art::InputTag fPMTWaveBaseLabel; + art::InputTag fPMTWaveSimLabel; + int fPMTWaveTestCh; + std::vector< art::InputTag > fPMTHitInputLabels; + + // CRT + bool fCRTOverlayHits; + std::vector< art::InputTag > fCRTHitInputLabels; + + // PROBABLY TEMPORARY + //geo::GeometryCore const* fGeom; +}; + + +OverlayProducts::OverlayProducts(fhicl::ParameterSet const& p) + : EDProducer{p}, + fTPCOverlayRaw ( p.get< bool >("TPCOverlayRaw", false) ), + fTPCOverlayROI ( p.get< bool >("TPCOverlayROI", false) ), + fTPCOverlayHits ( p.get< bool >("TPCOverlayHits", false) ), + fTPCHitsWireAssn ( p.get< bool >("TPCHitsWireAssn", true) ), + fTPCRawInputLabels ( p.get< std::vector >("TPCRawInputLabels", {}) ), + fTPCRawOutputLabel ( p.get< art::InputTag >("TPCRawOutputLabel","") ), + fTPCROIInputLabels ( p.get< std::vector >("TPCROIInputLabels", {}) ), + fTPCHitInputLabels ( p.get< std::vector >("TPCHitInputLabels", {}) ), + fTPCHitCreatorInstanceName ( p.get("TPCHitCreatorInstanaceName","") ), + fPMTOverlayRaw ( p.get< bool >("PMTOverlayRaw", false) ), + fPMTOverlayHits ( p.get< bool >("PMTOverlayHits", false) ), + fPMTWaveDataLabel ( p.get< art::InputTag >("PMTWaveDataLabel", "") ), + fPMTWaveBaseLabel ( p.get< art::InputTag >("PMTWaveBaseLabel", "") ), + fPMTWaveSimLabel ( p.get< art::InputTag >("PMTWaveSimLabel", "") ), + fPMTWaveTestCh ( p.get< int >("PMTWaveTestCh", -1) ), + fPMTHitInputLabels ( p.get< std::vector >("PMTHitInputLabels", {}) ), + fCRTOverlayHits ( p.get< bool >("CRTOverlayHits", false) ), + fCRTHitInputLabels ( p.get< std::vector >("CRTHitInputLabels", {}) ) + // More initializers here. +{ + // Call appropriate produces<>() functions here. + // Call appropriate consumes<>() for any products to be retrieved by this module. + + if ( !fTPCOverlayRaw && !fTPCOverlayROI && !fTPCOverlayHits && !fPMTOverlayRaw && !fPMTOverlayHits && !fCRTOverlayHits ) { + throw cet::exception("OverlayProducts") << "Error... need to be doing SOME overlay." << std::endl; + } + + if ( (fTPCOverlayRaw && fTPCRawInputLabels.size()<2) || + (fTPCOverlayROI && fTPCROIInputLabels.size()<2) || + (fTPCOverlayHits && fTPCHitInputLabels.size()<2) || + (fPMTOverlayRaw && fPMTWaveDataLabel==fPMTWaveSimLabel) || + (fPMTOverlayHits && fPMTHitInputLabels.size()<2) || + (fCRTOverlayHits && fCRTHitInputLabels.size()<2) ) + { + throw cet::exception("OverlayProducts") << "Error... Need at least 2 input products to merge." << std::endl; + } + + if ( fTPCOverlayRaw ) produces< std::vector >(fTPCRawOutputLabel.instance()); + if ( fTPCOverlayROI ) produces< std::vector >(); + // basically from the GausHitFinder + if ( fTPCOverlayHits ) recob::HitCollectionCreator::declare_products(producesCollector(), fTPCHitCreatorInstanceName, fTPCHitsWireAssn, false); + + if ( fPMTOverlayRaw ) { + produces< std::vector >(); + produces< std::vector >(); + produces< art::Assns >(); + } + if ( fPMTOverlayHits ) produces< std::vector >(); + + if ( fCRTOverlayHits ) produces< std::vector >(); +} + +void OverlayProducts::produce(art::Event& e) +{ + // Implementation of required member function here. + + // PROBABLY TEMPORARY + //fGeom = lar::providerFrom(); + + // TODO: Add CRT systems... + + // NOTE! IF OVERLAYING RAW DIGITS/WAVEFORMS, DO --DATA-- (OR JUST WHATEVER SAMPLE HAS THE NOISE AND THE BASE TIME) FIRST IN THE VECTOR OF LABELS... + + // TODO: recob::Hit mixing will just join the two products together for now. Probably need better version in future... + // NOTE: see also HitMerger_module.cc + + if ( fTPCOverlayRaw ) + { + // Need to see what Assns exist for the raw::RawDigits but this is a start... + std::map< raw::ChannelID_t, raw::RawDigit::ADCvector_t > rawdigitMap; + std::map< raw::ChannelID_t, float > pedestalMap; + std::map< raw::ChannelID_t, float > sigmaMap; + std::map< raw::ChannelID_t, raw::Compress_t > compressMap; + + for ( auto const& iLabel : fTPCRawInputLabels ) { + art::Handle< std::vector > digitsHandle; + std::vector< art::Ptr > digits; + if ( e.getByLabel(iLabel,digitsHandle) ) { + art::fill_ptr_vector(digits,digitsHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find raw::RawDigit with label " << iLabel << "."; + return; + } + + for ( auto const& iDigit : digits ) { + // Since the overlaid product is pure signal (no noise) we will keep the same pedestal and sigma, so only write to the map if no previous obect... + auto chID = iDigit->Channel(); + //std::vector wire = fGeom->ChannelToWire( chID ); + //std::cout << chID << " " << wire[0].Cryostat << " " << wire[0].TPC << " " << wire[0].Plane << " " << wire[0].Wire << std::endl; + if ( pedestalMap.find( chID ) == pedestalMap.end() ) { + pedestalMap[ chID ] = iDigit->GetPedestal(); + sigmaMap[ chID ] = iDigit->GetSigma(); + compressMap[ chID ] = iDigit->Compression(); + rawdigitMap[ chID ] = iDigit->ADCs(); + } + else { + auto const& thisADC = iDigit->ADCs(); + for ( unsigned int idx = 0; idx < rawdigitMap[ chID ].size(); ++idx ) { + rawdigitMap[ chID ][ idx ] += thisADC[ idx ]; + } + } + } // loop RawDigits + } // loop labels + + // Save the overlaid product + std::unique_ptr< std::vector< raw::RawDigit > > rawDigitVec = std::make_unique< std::vector< raw::RawDigit > >(); + + //unsigned int intCheck=0; + for ( auto const &[chKey, adcVals] : rawdigitMap ) { + const raw::ChannelID_t channel = chKey; + const raw::RawDigit::ADCvector_t ADCs = adcVals; + + const float pedestal = pedestalMap[ chKey ]; + const float sigma = sigmaMap[ chKey ]; + const raw::Compress_t compression = compressMap[ chKey ]; + + //std::cout << intCheck << " " << channel << " " << ADCs.size() << " " << compression << std::endl; + + raw::RawDigit rawDigit(channel, ADCs.size(), ADCs, compression); + rawDigit.SetPedestal( pedestal, sigma ); + rawDigitVec->push_back(rawDigit); + + //std::cout << " " << rawDigit.GetSigma() << std::endl; + //intCheck+=1; + } + e.put( std::move(rawDigitVec) ); + } + /* + if ( fTPCOverlayROI ){ + std::map< raw::ChannelID_t, std::vector< std::vector::size_type > > chOffsetsMap; + std::map< raw::ChannelID_t, std::vector< std::vector > > chDataMap; + + for ( auto const& iLabel : fTPCROIInputLabels ) { + art::Handle< std::vector > chROIsHandle; + std::vector< art::Ptr > chROIs; + if ( e.getByLabel(iLabel,chROIsHandle) ) { + art::fill_ptr_vector(chROIs,chROIsHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find recob::ChannelROIs with label " << iLabel << "."; + return; + } + + for ( auto const& ichROI : chROIs ) { + //std::cout << "New ChannelROI" << std::endl; + const recob::ChannelROI::RegionsOfInterest_t& channelROIs = ichROI->SignalROI(); + const raw::ChannelID_t chid = ichROI->Channel(); + + for ( auto const& range : channelROIs.get_ranges() ) { + //size_t startTick = range.begin_index(); + //std::vector dataVec(range.data().size()); + //for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) dataVec[binIdx] = range.data()[binIdx]; + //std::cout << " range, start: " << startTick << ", " << dataVec.size() << " size -- [0] = " << dataVec[0] << std::endl; + + std::vector dataVec(range.data().size()); + for(size_t binIdx=0; binIdx < range.data().size(); ++binIdx) dataVec[binIdx] = range.data()[binIdx]; + + // If channel not in the map already then add it to the map. Otherwise, add the ROI to the map, checking if it overlaps with an already present one... + if ( chOffsetsMap.find( chid ) == chOffsetsMap.end() ) { + chOffsetsMap[ chid ] = std::vector< std::vector::size_type >{ range.begin_index() }; + chDataMap[ chid ] = std::vector< std::vector >{ dataVec }; + } + else { + // Check if this ROI overlaps with any other region for this channel + std::vector::size_type loRange = range.begin_index(); + std::vector::size_type hiRange = range.begin_index() + dataVec.size(); + + std::vector< unsigned int > overlaps; + + for( unsigned int iOffset=0; iOffset < chOffsetsMap[chid].size(); ++iOffset ) { + if ( ( loRange > chOffsetsMap[chid].at(iOffset) && loRange < chOffsetsMap[chid].at(iOffset)+chDataMap[chid].at(iOffset).size() ) || + ( hiRange > chOffsetsMap[chid].at(iOffset) && hiRange < chOffsetsMap[chid].at(iOffset)+chDataMap[chid].at(iOffset).size() ) ) { + overlaps.push_back( iOffset ); + } + } + + // If no overlaps, put the ROI into the map and move on + if ( overlaps.size() == 0 ) { + chOffsetsMap[ chid ].push_back( loRange ); + chDataMap[ chid ].push_back( dataVec ); + continue; + } + + // If overlaps, find the min offset and max offset+size to appropriately pad out the vectors + unsigned int idxMinOffset = 0; + std::vector::size_type minOffset = loRange; + std::vector::size_type maxEnd = hiRange; + for ( auto const& overlap : overlaps ) { + if ( chOffsetsMap[chid].at(overlap) < minOffset ){ + minOffset = chOffsetsMap[chid].at(overlap); + idxMinOffset = overlap; + } + if ( maxEnd < chOffsetsMap[chid].at(overlap)+chDataMap[chid].at(overlap).size() ) maxEnd = chOffsetsMap[chid].at(overlap)+chDataMap[chid].at(overlap).size(); + } + + // Sum up the ROIs + // 1. pad this data vector + for ( unsigned int iBegin=0; iBegin < loRange-minOffset; ++iBegin) dataVec.insert( dataVec.begin(), 0 ); + for ( unsigned int iEnd=0; iEnd < maxEnd-hiRange; ++iEnd) dataVec.push_back( 0 ); + // 2. now add in the other vectors + for ( auto const& overlap : overlaps ) { + std::vector::size_type delta = chOffsetsMap[chid].at(overlap) - minOffset; + for ( unsigned int idxData=0; idxData > channelROIVec = std::make_unique< std::vector< recob::ChannelROI > >(); + + for ( auto const &[chKey, offsetVec] : chOffsetsMap ) { + recob::ChannelROI::RegionsOfInterest_t ROIVec; + for ( unsigned int idxROI = 0; idxROI < offsetVec.size(); ++idxROI ) { + ROIVec.add_range(offsetVec[idxROI], chDataMap[chKey].at(idxROI) ); + } + + recob::ChannelROI thisChROI(ROIVec,chKey); + channelROIVec->push_back(thisChROI); + } + + e.put( std::move(channelROIVec) ); + }*/ + + if ( fTPCOverlayHits ) + { + // As in Gauss Hit Finder code (larreco/HitFinder/GausHitFinder_module.cc) + recob::HitCollectionCreator hitCol(e, fTPCHitCreatorInstanceName, fTPCHitsWireAssn, false); + + // Loop through hit labels and put together the hit set + for ( auto const& iLabel : fTPCHitInputLabels ) { + art::Handle< std::vector > hitsHandle; + std::vector< art::Ptr > hits; + if ( e.getByLabel(iLabel,hitsHandle) ) { + art::fill_ptr_vector(hits,hitsHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find recob::Hit with label " << iLabel << "."; + return; + } + + art::FindManyP fmwire(hitsHandle, e, iLabel); + if( !fmwire.isValid() && fTPCHitsWireAssn ){ + mf::LogError("OverlayProducts") << "Error in validity of fmwire. Returning."; + return; + } + + for ( auto const& iHitPtr : hits ) { + recob::Hit theHit = *iHitPtr; + // And get the associated wire -- if we should + if ( fTPCHitsWireAssn ) { + std::vector< art::Ptr > hitWires = fmwire.at(iHitPtr.key()); + if ( hitWires.size() == 0 ) throw cet::exception("OverlayProducts") << "Hit found with no associated wires...\n"; + else if ( hitWires.size() > 1 ) mf::LogWarning("OverlayProducts") << "Hit with >1 recob::Wire associated..."; + hitCol.emplace_back(theHit,hitWires[0]); + } + else hitCol.emplace_back(theHit); + } + } + + // Put the hit into the event + hitCol.put_into(e); + } + + // PMT Overlays + // Get the data waveforms and populate the map. Then get the sim waveforms and put them into the right waveforms. + if ( fPMTOverlayRaw ) { + // Produced vectors + std::unique_ptr< std::vector< raw::OpDetWaveform > > opWaveformVec = std::make_unique< std::vector< raw::OpDetWaveform > >(); + std::unique_ptr< std::vector< icarus::WaveformBaseline > > baselineVec = std::make_unique< std::vector< icarus::WaveformBaseline > >(); + art::Assns baselineToWaveforms; + + std::map< raw::Channel_t, std::vector< opWaveformOverlay > > opWaveformsMap; + // Load in data waveforms + art::Handle< std::vector > dataWavesHandle; + std::vector< art::Ptr > dataWaves; + if ( e.getByLabel(fPMTWaveDataLabel,dataWavesHandle) ) { + art::fill_ptr_vector(dataWaves,dataWavesHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find raw::OpDetWaveform with label " << fPMTWaveDataLabel << "."; + return; + } + // ... and the baselines + art::Handle< std::vector > baselinesHandle; + std::vector< art::Ptr > baselines; + if ( e.getByLabel(fPMTWaveBaseLabel,baselinesHandle) ) { + art::fill_ptr_vector(baselines,baselinesHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find raw::OpDetWaveform with label " << fPMTWaveBaseLabel << "."; + return; + } + art::FindManyP fmbase( dataWavesHandle, e, fPMTWaveBaseLabel ); + if( !fmbase.isValid() ) { + mf::LogError("OverlayProducts") << "Error in validity of fmbase. Returning."; + return; + } + // TEST!!!! PRINT BASELINES //////////////// + std::cout << "// ---- BASELINES ---- " << std::endl; + for ( unsigned int idxBase=0; idxBasebaseline() << " "; + } + std::cout << "---- BASELINES ---- //" << std::endl; + std::cout << std::endl; + //////////////////////////////////////////// + // Now fill the map + for ( auto const& dataWave : dataWaves ) { + auto chID = dataWave->ChannelNumber(); + auto time = dataWave->TimeStamp(); + auto wvfm = dataWave->Waveform(); + + if ( wvfm.size() == 0 ) continue; + + if ( opWaveformsMap.find( chID ) == opWaveformsMap.end() ) + opWaveformsMap[ chID ] = std::vector< opWaveformOverlay >(); + + auto baselines = fmbase.at( dataWave.key() ); + if ( baselines.size() > 1 ) std::cout << "SHOULD NOT BE MORE THAN ONE BASELINE..." << std::endl; + + opWaveformOverlay dataOpWave; + dataOpWave.timestamp = time; + dataOpWave.wvfm = wvfm; + dataOpWave.baseline = ( baselines.size()>0 ? baselines[0]->baseline() : wvfm[0] ); // TODO: better way of guarding? + + // Print + if ( int(chID)==fPMTWaveTestCh && wvfm.size()>10500 ) { + std::cout << "DATA WAVEFORM AT TIME " << std::setprecision(8) << time << ": [ "; + for ( unsigned int i=0; i > simWavesHandle; + std::vector< art::Ptr > simWaves; + if ( e.getByLabel(fPMTWaveSimLabel,simWavesHandle) ) { + art::fill_ptr_vector(simWaves,simWavesHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find raw::OpDetWaveform with label " << fPMTWaveSimLabel << "."; + return; + } + for ( auto const& simWave : simWaves ) { + auto chID = simWave->ChannelNumber(); + auto time = simWave->TimeStamp(); + auto wvfm = simWave->Waveform(); + raw::ADC_Count_t simBaseline = 14999; // TODO: automate, parameter, not hard-code? + + if ( wvfm.size() == 0 ) continue; + + // PRINT + if ( int(chID)==fPMTWaveTestCh && wvfm.size()>10500 ) { + std::cout << "SIM WAVEFORM AT TIME " << std::setprecision(8) << time << ": [ "; + for ( unsigned int i=0;i overlaps; + if ( int(chID)==fPMTWaveTestCh && wvfm.size()>10500 ) { + std::cout << "Overlap idx vals: [ "; + } + for ( unsigned int idxWave=0; idxWave thisWaveformOverlay.timestamp ) || + ( time < thisWaveformOverlay.endTime() && simEnd > thisWaveformOverlay.endTime() ) ) { + overlaps.push_back( idxWave ); + } + + if ( int(chID)==fPMTWaveTestCh && wvfm.size()>10500 ) { + std::cout << thisWaveformOverlay.timestamp << "-to-" << thisWaveformOverlay.endTime() << " "; + } + } // iterate saved waveforms looking for overlap + if ( int(chID)==fPMTWaveTestCh && wvfm.size()>10500 ) { + std::cout << " ]" << std::endl; + } + + if ( int(chID)==fPMTWaveTestCh && wvfm.size()>10500 ) { + std::cout << std::endl; + std::cout << "OVERLAPS = " << overlaps.size() << std::endl; + std::cout << " [ "; + for ( auto const& i : overlaps ) std::cout << i << " "; + std::cout << "]" << std::endl; + std::cout << std::endl; + } + + // Case 1: No overlaps. Make a new entry in the map as above... + if ( overlaps.size() == 0 ) { + opWaveformOverlay simOpWave; + simOpWave.timestamp = time; + simOpWave.wvfm = wvfm; + simOpWave.baseline = simBaseline; + opWaveformsMap[chID].push_back(simOpWave); + } + // Case 2: 1+ overlaps. Make a new struct and erase the overlaps + if ( overlaps.size() >= 1 ) { + // Make a map from timestamp to idx within overlap to find the minimum time and then break these out into vectors + // TODO: better way? + std::map< raw::TimeStamp_t, unsigned int > timestampIdxMap; + for ( unsigned int idxOverlap=0; idxOverlap orderedTimes; + std::vector< unsigned int > orderedOverlaps; + for ( auto const &[time, idx] : timestampIdxMap ) { + orderedTimes.push_back(time); + orderedOverlaps.push_back(idx); + } + + if ( int(chID)==fPMTWaveTestCh && wvfm.size()>10500 ) { + std::cout << "Sim Time: " << time << std::endl; + std::cout << "Times: [ "; + for ( unsigned int i=0; i adcVec; + double deltaT = opWaveformsMap[ chID ].at( orderedOverlaps[0] ).deltaT; + + unsigned int idxWvfmEntry=0; + double loTime=std::numeric_limits::lowest(); + for ( unsigned int iOverlap=0; iOverlap loTime && time+((double(idxWvfmEntry)+0.5)*deltaT) < orderedTimes[iOverlap] ) { + adcVec.push_back( overlayOpWave.baseline + (wvfm[idxWvfmEntry] - simBaseline) ); + idxWvfmEntry+=1; + } + + opWaveformOverlay thisOverlap = opWaveformsMap[ chID ].at( orderedOverlaps[iOverlap] ); + + // If after, then start putting in entries before the overlapping time + unsigned int alreadyInOverlap = 0; + while ( idxWvfmEntry10500 ) { + std::cout << "... We should now erase entry " << positionToErase << " which is STRUCT WITH:" << std::endl; + std::cout << " time: " << opWaveformsMap[ chID ].at( positionToErase+1 ).timestamp << std::endl; + std::cout << " baseline: " << opWaveformsMap[ chID ].at( positionToErase+1 ).baseline << std::endl; + std::cout << " wv size: " << opWaveformsMap[ chID ].at( positionToErase+1 ).wvfm.size() << std::endl; + std::cout << " end time: " << opWaveformsMap[ chID ].at( positionToErase+1 ).endTime() << std::endl; + } + + opWaveformsMap[ chID ].erase( opWaveformsMap[ chID ].begin()+positionToErase ); + } + + // And save this new one + overlayOpWave.wvfm = adcVec; + opWaveformsMap[ chID ].push_back( overlayOpWave ); + } // deal with overlaps + } // look for overlaps + } // loop simulated waveforms + + // Now put all the waveforms into the collections + art::PtrMaker opdetwvfmPtrMaker(e); + art::PtrMaker baselinePtrMaker(e); + for ( auto const &[chID, wvfmStructVec] : opWaveformsMap ) { + for ( auto const& wvfmStruct : wvfmStructVec ) { + const raw::TimeStamp_t time = wvfmStruct.timestamp; + const raw::Channel_t chan = chID; + std::vector< raw::ADC_Count_t > vals = wvfmStruct.wvfm; + const icarus::WaveformBaseline::Baseline_t base = wvfmStruct.baseline; + + // Print + if ( int(chID)==fPMTWaveTestCh ) { + std::cout << "Saving STRUCT WITH: " << std::endl; + std::cout << " time: " << wvfmStruct.timestamp << std::endl; + std::cout << " baseline: " << wvfmStruct.baseline << std::endl; + std::cout << " wv size: " << wvfmStruct.wvfm.size() << std::endl; + std::cout << " end time: " << wvfmStruct.endTime() << std::endl; + } + + if ( int(chID)==fPMTWaveTestCh && vals.size()>10500 ) { + std::cout << "OVERLAID WAVEFORM AT TIME " << std::setprecision(8) << time << ": [ "; + for ( unsigned int i=0;i valsUnsigned; + for ( auto const& val : vals ) valsUnsigned.push_back( val ); + + raw::OpDetWaveform odw( time, chan, valsUnsigned ); + icarus::WaveformBaseline bLine( base ); + + opWaveformVec->push_back( odw ); + baselineVec->push_back( bLine ); + baselineToWaveforms.addSingle( baselinePtrMaker(baselineVec->size()-1), opdetwvfmPtrMaker(opWaveformVec->size()-1) ); + } + } + + // ... and put collections into the event + e.put( std::move( opWaveformVec ) ); + e.put( std::move( baselineVec ) ); + e.put( std::move( std::make_unique< art::Assns >(std::move(baselineToWaveforms)) ) ); + } + + if ( fPMTOverlayHits ) { + std::unique_ptr< std::vector< recob::OpHit > > opHitVec = std::make_unique< std::vector< recob::OpHit > >(); + + for ( auto const& iLabel : fPMTHitInputLabels ) { + art::Handle< std::vector > hitsHandle; + std::vector< art::Ptr > hits; + if ( e.getByLabel(iLabel,hitsHandle) ) { + art::fill_ptr_vector(hits,hitsHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find recob::OpHit with label " << iLabel << "."; + return; + } + + for ( auto const& iHit : hits ) { + recob::OpHit newOpHit( iHit->OpChannel(), iHit->PeakTime(), iHit->PeakTimeAbs(), iHit->Frame(), iHit->Width(), iHit->Area(), iHit->Amplitude(), iHit->PE(), iHit->FastToTotal() ); + opHitVec->push_back( newOpHit ); + } // hit + } // label + + e.put( std::move(opHitVec) ); + } + + // CRT Overlays + if ( fCRTOverlayHits ) { + std::unique_ptr< std::vector< sbn::crt::CRTHit > > crtHitVec = std::make_unique< std::vector< sbn::crt::CRTHit > >(); + + for ( auto const& iLabel : fCRTHitInputLabels ) { + art::Handle< std::vector > hitsHandle; + std::vector< art::Ptr > hits; + if ( e.getByLabel(iLabel,hitsHandle) ) { + art::fill_ptr_vector(hits,hitsHandle); + } + else{ + mf::LogWarning("OverlayProducts") << "Event failed to find sbn::crt::CRTHit with label " << iLabel << "."; + return; + } + + for ( auto const& iHit : hits ) { + sbn::crt::CRTHit newCRTHit; + newCRTHit.feb_id = iHit->feb_id; + newCRTHit.pesmap = iHit->pesmap; + newCRTHit.peshit = iHit->peshit; + newCRTHit.ts0_s = iHit->ts0_s; + newCRTHit.ts0_s_corr = iHit->ts0_s_corr; + newCRTHit.ts0_ns = iHit->ts0_ns; + newCRTHit.ts0_ns_corr = iHit->ts0_ns_corr; + newCRTHit.ts1_ns = iHit->ts1_ns; + newCRTHit.plane = iHit->plane; + newCRTHit.x_pos = iHit->x_pos; + newCRTHit.x_err = iHit->x_err; + newCRTHit.y_pos = iHit->y_pos; + newCRTHit.y_err = iHit->y_err; + newCRTHit.z_pos = iHit->z_pos; + newCRTHit.z_err = iHit->z_err; + newCRTHit.tagger = iHit->tagger; + + crtHitVec->push_back( newCRTHit ); + } // hit + } // label + + e.put( std::move(crtHitVec) ); + } +} + +DEFINE_ART_MODULE(OverlayProducts) diff --git a/icaruscode/Overlays/CMakeLists.txt b/icaruscode/Overlays/CMakeLists.txt new file mode 100644 index 000000000..63088399c --- /dev/null +++ b/icaruscode/Overlays/CMakeLists.txt @@ -0,0 +1,52 @@ +find_package(gallery REQUIRED) + +#find_package(larwirecell REQUIRED ) +#find_package(jsoncpp REQUIRED) +#find_package(spdlog REQUIRED) + +include_directories( $ENV{JSONCPP_INC} ) +link_directories( $ENV{JSONCPP_LIB} ) +include_directories( $ENV{SPDLOG_INC} ) +link_directories( $ENV{SPDLOG_LIB} ) +include_directories( $ENV{WIRECELL_INC} ) +link_directories( $ENV{WIRECELL_LIB} ) +include_directories( $ENV{LARWIRECELL_INC} ) +link_directories( $ENV{LARWIRECELL_LIB} ) + +cet_build_plugin( + SimInfoOverlayFilter art::EDFilter + LIBRARIES + PRIVATE + larcoreobj::SummaryData + lardataobj::Simulation + nusimdata::SimulationBase + gallery::gallery +) + +# This is what's needed for building the WireCellubsim WCT component +art_make_library( + LIBRARY_NAME WireCellICARUSDrifter + SOURCE ICARUSDrifter.cxx + LIBRARIES + art::Framework_Services_Registry + art::Utilities + ${Boost_SYSTEM_LIBRARY} + ${ROOT_BASIC_LIB_LIST} + ${JSONCPP} + larevt::CalibrationDBI_Providers + lardata::ArtDataHelper + ${FHICLCPP} + cetlib_except::cetlib_except + WireCellGen + WireCellUtil + WireCellAux + WireCellIface + WireCellPgraph + WireCellRoot + jsoncpp +) + + +install_headers() +install_source() +install_fhicl() diff --git a/icaruscode/Overlays/ICARUSDrifter.cxx b/icaruscode/Overlays/ICARUSDrifter.cxx new file mode 100644 index 000000000..741e57036 --- /dev/null +++ b/icaruscode/Overlays/ICARUSDrifter.cxx @@ -0,0 +1,77 @@ +#include "ICARUSDrifter.h" +#include "art/Framework/Principal/Event.h" +#include "larevt/CalibrationDBI/Providers/DBFolder.h" +#include "larcore/Geometry/Geometry.h" +#include "WireCellUtil/NamedFactory.h" +#include "WireCellUtil/Units.h" +#include "WireCellUtil/Persist.h" +#include + +WIRECELL_FACTORY(wclsICARUSDrifter, wcls::ICARUSDrifter, + wcls::IArtEventVisitor, WireCell::IDrifter) + +using namespace WireCell; + +wcls::ICARUSDrifter::ICARUSDrifter() + : Drifter(), + fDB(NULL) +{ +} + +wcls::ICARUSDrifter::~ICARUSDrifter() +{ + if (fDB) delete fDB; +} + +double wcls::ICARUSDrifter::GetLifetime(uint64_t run, uint64_t itpc) { + // check the cache + if (fLifetimes.count({run, itpc})) { + if (fVerbose) std::cout << "Run: " << run << " ITPC: " << itpc << " SEEN. Returning: " << fLifetimes.at({run, itpc}) << ".\n"; + return fLifetimes.at({run, itpc}); + } + + // Look up the run + // + // Translate the run into a fake "timestamp" + fDB->UpdateData((run+1000000000)*1000000000); + + double tau; + fDB->GetNamedChannelData(itpc, "elifetime", tau); + + // Set the cache + std::pair toset {run, itpc}; + fLifetimes[toset] = tau; + + if (fVerbose) std::cout << "Run: " << run << " ITPC: " << itpc << " NEW. Returning: " << tau << ".\n"; + + return tau; +} + +void wcls::ICARUSDrifter::visit(art::Event & event) +{ + /// this function will be executed post WCT configuration !!! + if (fVerbose) std::cout << "wcls ICARUSDrifter: electron lifetime read in!\n"; + + if (fELifetimeCorrection) { + double elifetime = GetLifetime(event.id().runID().run(), fTPC); + Drifter::set_lifetime(elifetime*units::us); + + if (fVerbose) std::cout << "www: lifetime from database = " << elifetime << std::endl; + if (fVerbose) std::cout << "www: lifetime from database = " << elifetime*units::us << std::endl; + } + else { + Drifter::set_lifetime(1000.0*units::ms); // set to infinite + } +} + +void wcls::ICARUSDrifter::configure(const WireCell::Configuration& cfg) +{ + Drifter::configure(cfg); + + fDBFileName = cfg["DBFileName"].asString(); + fDBTag = cfg["DBTag"].asString(); + fVerbose = cfg["Verbose"].asBool(); + fDB = new lariov::DBFolder(fDBFileName, "", "", fDBTag, true, false); + fTPC = cfg["TPC"].asInt(); + fELifetimeCorrection = cfg["ELifetimeCorrection"].asBool(); +} diff --git a/icaruscode/Overlays/ICARUSDrifter.h b/icaruscode/Overlays/ICARUSDrifter.h new file mode 100644 index 000000000..a490eb1a3 --- /dev/null +++ b/icaruscode/Overlays/ICARUSDrifter.h @@ -0,0 +1,49 @@ +#ifndef ICARUS_DRIFTER_H +#define ICARUS_DRIFTER_H + +#include "WireCellGen/Drifter.h" + +#include "larwirecell/Interfaces/IArtEventVisitor.h" +#include "larevt/CalibrationDBI/Providers/DBFolder.h" + +#include +#include + +namespace wcls { + class ICARUSDrifter : public WireCell::Gen::Drifter, + public IArtEventVisitor + { + public: + ICARUSDrifter(); + virtual ~ICARUSDrifter(); + + /// IArtEventVisitor. + // + // Note: we don't actually poke at the event but use this + // entry to refresh info from services in case they change + // event-to-event. + virtual void visit(art::Event & event); + + /// IConfigurable. + // + // Defer default to parent. By default this class does not + // overriding. + + // Defer to parent but override if asked to. + virtual void configure(const WireCell::Configuration& config); + + private: + double GetLifetime(uint64_t run, uint64_t itpc); + + std::string fDBFileName; + std::string fDBTag; + bool fVerbose; + lariov::DBFolder *fDB; + + int fTPC; + bool fELifetimeCorrection; + + std::map, double> fLifetimes; // cache + }; +} // wcls +#endif diff --git a/icaruscode/Overlays/SimInfoOverlayFilter_module.cc b/icaruscode/Overlays/SimInfoOverlayFilter_module.cc new file mode 100644 index 000000000..e4fec2ab3 --- /dev/null +++ b/icaruscode/Overlays/SimInfoOverlayFilter_module.cc @@ -0,0 +1,803 @@ +//////////////////////////////////////////////////////////////////////// +// Class: SimInfoOverlayFilter +// Plugin Type: filter (art v2_11_03) +// File: SimInfoOverlayFilter_module.cc +// +// Generated at Sat Nov 24 16:09:21 2018 by Wesley Ketchum using cetskelgen +// from cetlib version v3_03_01. +//////////////////////////////////////////////////////////////////////// + +#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/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "canvas/Persistency/Provenance/ProductID.h" +#include "canvas/Persistency/Provenance/canonicalProductName.h" +#include "canvas/Utilities/FriendlyName.h" +#include "art/Persistency/Common/PtrMaker.h" + +#include +#include +#include +#include +#include +#include + +//gallery... +#include "gallery/Event.h" +#include "gallery/Handle.h" +#include "gallery/ValidHandle.h" + +//association includes +#include "canvas/Persistency/Common/FindMany.h" +#include "canvas/Persistency/Common/FindOne.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "canvas/Persistency/Common/Assns.h" + +//generator products +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/GTruth.h" +#include "nusimdata/SimulationBase/MCFlux.h" +#include "lardataobj/Simulation/BeamGateInfo.h" + +//G4 products +#include "nusimdata/SimulationBase/MCParticle.h" +#include "lardataobj/Simulation/SimPhotons.h" +#include "lardataobj/Simulation/SimChannel.h" +#include "lardataobj/Simulation/AuxDetSimChannel.h" +#include "lardataobj/Simulation/SimEnergyDeposit.h" +#include "lardataobj/Simulation/GeneratedParticleInfo.h" +#include "lardataobj/MCBase/MCParticleLite.h" + +//pot product +#include "larcoreobj/SummaryData/POTSummary.h" + + +namespace mix { + class SimInfoOverlayFilter; +} + + +class mix::SimInfoOverlayFilter : public art::EDFilter { +public: + explicit SimInfoOverlayFilter(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. + SimInfoOverlayFilter(SimInfoOverlayFilter const &) = delete; + SimInfoOverlayFilter(SimInfoOverlayFilter &&) = delete; + SimInfoOverlayFilter & operator = (SimInfoOverlayFilter const &) = delete; + SimInfoOverlayFilter & operator = (SimInfoOverlayFilter &&) = delete; + + // Required functions. + bool filter(art::Event & e) override; + + // Selected optional functions. + void beginJob() override; + bool beginRun(art::Run & r) override; + bool beginSubRun(art::SubRun & sr) override; + void endJob() override; + bool endRun(art::Run & r) override; + bool endSubRun(art::SubRun & sr) override; + void respondToCloseInputFile(art::FileBlock const & fb) override; + void respondToCloseOutputFiles(art::FileBlock const & fb) override; + void respondToOpenInputFile(art::FileBlock const &fb) override; + void respondToOpenOutputFiles(art::FileBlock const & fb) override; + +private: + + //Infrastructure to read the event... + std::vector fSimInputFileNames; + gallery::Event gEvent; + + //input tag lists... + //note, if these are empty, that's ok: then we won't produce them on the event + + std::vector fMCTruthInputModuleLabels; + std::vector fMCFluxInputModuleLabels; + std::vector fGTruthInputModuleLabels; + std::vector fBeamGateInputModuleLabels; + + std::vector fMCParticleInputModuleLabels; + std::vector fSimEnergyDepositInputModuleLabels; + std::vector fAuxDetSimChannelInputModuleLabels; + std::vector fMCParticleLiteInputModuleLabels; + std::vector fSimChannelInputModuleLabels; + std::vector fSimPhotonsInputModuleLabels; + + //now for the associations... + std::vector fMCTruthMCFluxAssnsInputModuleLabels; + std::vector fMCTruthGTruthAssnsInputModuleLabels; + std::vector fMCTruthMCParticleAssnsInputModuleLabels; + + // Added by C Thorpe: New fhicl parameters needed to handle separate gen/g4 stages in auxiallary file + bool fAuxHasSeparateGenG4; + std::vector fMCTruthMCParticleAssnsMCTruthLookupLabel; + + template + using CollectionMap = std::unordered_map< std::string, std::unique_ptr >; + + CollectionMap< std::vector > fMCTruthMap; + CollectionMap< std::vector > fMCFluxMap; + CollectionMap< std::vector > fGTruthMap; + CollectionMap< std::vector > fBeamGateInfoMap; + + CollectionMap< art::Assns > fMCTruthMCFluxAssnsMap; + CollectionMap< art::Assns > fMCTruthGTruthAssnsMap; + + CollectionMap< std::vector > fMCParticleMap; + + CollectionMap< std::vector > fSimEnergyDepositMap; + CollectionMap< std::vector > fAuxDetSimChannelMap; + CollectionMap< std::vector > fMCParticleLiteMap;//Ivan + CollectionMap< std::vector > fSimChannelMap; + CollectionMap< std::vector > fSimPhotonsMap; + + CollectionMap< art::Assns > fMCTruthMCParticleAssnsMap; + //verbose + int fVerbosity; + + //pot accounting + bool fFillPOTInfo; + art::InputTag fPOTSummaryTag; + double fPOTSum_totpot; + double fPOTSum_totgoodpot; + double fPOTSum_totspills; + double fPOTSum_goodspills; + + void MakePOTMap(); + std::map fSR_POTPerEvent; + std::map fSR_GoodPOTPerEvent; + std::map fSR_SpillsPerEvent; + std::map fSR_GoodSpillsPerEvent; + + void FillInputModuleLabels(fhicl::ParameterSet const &, + std::string, + std::vector &); + void FillInputModuleLabels(fhicl::ParameterSet const &); + + template + void DeclareProduces(std::vector const&); + void DeclareProduces(); + + template + void InitializeCollectionMap(std::vector const&, + CollectionMap &); + void InitializeCollectionMaps(); + + template + void PutCollectionOntoEvent(art::Event &, + CollectionMap &); + void PutCollectionsOntoEvent(art::Event &); + + template + std::map< std::pair , art::Ptr > + FillCollectionMap(std::vector const&, + CollectionMap< std::vector > &, + std::string, + art::Event &); + + template< class T, class U > + void FillAssnsCollectionMap(std::vector const&, + CollectionMap< art::Assns > &, + std::map< std::pair , art::Ptr > const&, + std::map< std::pair , art::Ptr > const&); + template< class T, class U, class D> + void FillAssnsCollectionMap(std::vector const&, + CollectionMap< art::Assns > &, + std::map< std::pair , art::Ptr > const&, + std::map< std::pair , art::Ptr > const&); +}; + + + +mix::SimInfoOverlayFilter::SimInfoOverlayFilter(fhicl::ParameterSet const & p) + : EDFilter{p}, + fSimInputFileNames(p.get>("SimInputFileNames")), + gEvent(fSimInputFileNames), + fAuxHasSeparateGenG4(p.get("AuxHasSeparateGenG4",false)), + fVerbosity(p.get("Verbosity",-1)) +{ + FillInputModuleLabels(p); + DeclareProduces(); + + fFillPOTInfo = p.get("FillPOTInfo",true); + if(fFillPOTInfo){ + fPOTSummaryTag = p.get("POTSummaryTag"); + fPOTSum_totpot = 0.0; + fPOTSum_totgoodpot = 0.0; + fPOTSum_totspills = 0.0; + fPOTSum_goodspills = 0.0; + + this->produces(); + MakePOTMap(); + } + +} + +void mix::SimInfoOverlayFilter::FillInputModuleLabels(fhicl::ParameterSet const & p, + std::string s, + std::vector & labels) +{ + labels = p.get< std::vector >(s,std::vector()); +} + +void mix::SimInfoOverlayFilter::FillInputModuleLabels(fhicl::ParameterSet const & p) +{ + FillInputModuleLabels(p,"MCTruthInputModuleLabels",fMCTruthInputModuleLabels); + FillInputModuleLabels(p,"MCFluxInputModuleLabels",fMCFluxInputModuleLabels); + FillInputModuleLabels(p,"GTruthInputModuleLabels",fGTruthInputModuleLabels); + FillInputModuleLabels(p,"BeamGateInputModuleLabels",fBeamGateInputModuleLabels); + + FillInputModuleLabels(p,"MCParticleInputModuleLabels",fMCParticleInputModuleLabels); + FillInputModuleLabels(p,"SimEnergyDepositInputModuleLabels",fSimEnergyDepositInputModuleLabels); + FillInputModuleLabels(p,"AuxDetSimChannelInputModuleLabels",fAuxDetSimChannelInputModuleLabels); + FillInputModuleLabels(p,"MCParticleLiteInputModuleLabels",fMCParticleLiteInputModuleLabels); //Ivan + FillInputModuleLabels(p,"SimChannelInputModuleLabels",fSimChannelInputModuleLabels); + FillInputModuleLabels(p,"SimPhotonsInputModuleLabels",fSimPhotonsInputModuleLabels); + + FillInputModuleLabels(p,"MCTruthMCFluxAssnsInputModuleLabels",fMCTruthMCFluxAssnsInputModuleLabels); + FillInputModuleLabels(p,"MCTruthGTruthAssnsInputModuleLabels",fMCTruthGTruthAssnsInputModuleLabels); + FillInputModuleLabels(p,"MCTruthMCParticleAssnsInputModuleLabels",fMCTruthMCParticleAssnsInputModuleLabels); + + + // New fhicl parameter needed to handle separate gen/g4 passes in the aux file, crashes when trying to construct + // assns otherwise + + if(fAuxHasSeparateGenG4){ + FillInputModuleLabels(p,"MCTruthMCParticleAssnsMCTruthLookupLabel",fMCTruthMCParticleAssnsMCTruthLookupLabel); + if(!fMCTruthMCParticleAssnsMCTruthLookupLabel.size()) + throw cet::exception("SimInfoOverlayFilter") + << "Expecting aux file produced with separate gen/g4 passes, MCTruthMCParticleAssnsMCTruthLookupLabel in config needs to point to gen info in aux file"; + } +} + +template +void mix::SimInfoOverlayFilter::DeclareProduces(std::vector const& labels) +{ + for(auto label : labels) + this->produces(label.instance()); +} + +void mix::SimInfoOverlayFilter::DeclareProduces() +{ + DeclareProduces< std::vector >(fMCTruthInputModuleLabels); + DeclareProduces< std::vector >(fMCFluxInputModuleLabels); + DeclareProduces< std::vector >(fGTruthInputModuleLabels); + DeclareProduces< std::vector >(fBeamGateInputModuleLabels); + + DeclareProduces< std::vector >(fMCParticleInputModuleLabels); + DeclareProduces< std::vector >(fSimEnergyDepositInputModuleLabels); + DeclareProduces< std::vector >(fAuxDetSimChannelInputModuleLabels); + DeclareProduces< std::vector >(fMCParticleLiteInputModuleLabels); //Ivan + DeclareProduces< std::vector >(fSimChannelInputModuleLabels); + DeclareProduces< std::vector >(fSimPhotonsInputModuleLabels); + + DeclareProduces< art::Assns > + (fMCTruthMCFluxAssnsInputModuleLabels); + DeclareProduces< art::Assns > + (fMCTruthGTruthAssnsInputModuleLabels); + DeclareProduces< art::Assns > + (fMCTruthMCParticleAssnsInputModuleLabels); +} + +template +void mix::SimInfoOverlayFilter::InitializeCollectionMap(std::vector const& labels, + CollectionMap & colmap) +{ + for(auto const& label : labels) + colmap[label.instance()] = std::make_unique(); +} + +void mix::SimInfoOverlayFilter::InitializeCollectionMaps() +{ + InitializeCollectionMap(fMCTruthInputModuleLabels,fMCTruthMap); + InitializeCollectionMap(fMCFluxInputModuleLabels,fMCFluxMap); + InitializeCollectionMap(fGTruthInputModuleLabels,fGTruthMap); + InitializeCollectionMap(fBeamGateInputModuleLabels,fBeamGateInfoMap); + + InitializeCollectionMap(fMCTruthMCFluxAssnsInputModuleLabels,fMCTruthMCFluxAssnsMap); + InitializeCollectionMap(fMCTruthGTruthAssnsInputModuleLabels,fMCTruthGTruthAssnsMap); + + InitializeCollectionMap(fMCParticleInputModuleLabels,fMCParticleMap); + InitializeCollectionMap(fSimEnergyDepositInputModuleLabels,fSimEnergyDepositMap); + InitializeCollectionMap(fAuxDetSimChannelInputModuleLabels,fAuxDetSimChannelMap); + InitializeCollectionMap(fMCParticleLiteInputModuleLabels,fMCParticleLiteMap); //Ivan + InitializeCollectionMap(fSimChannelInputModuleLabels,fSimChannelMap); + InitializeCollectionMap(fSimPhotonsInputModuleLabels,fSimPhotonsMap); + + std::cout << "Filling fMCTruthMCParticleAssnsMap" << std::endl; + std::cout << "fMCTruthMCParticleAssnsInputModuleLabels.size()=" << fMCTruthMCParticleAssnsInputModuleLabels.size() << std::endl; + if(fMCTruthMCParticleAssnsInputModuleLabels.size()){ + std::cout << "Using label " << fMCTruthMCParticleAssnsInputModuleLabels.at(0) << std::endl; + auto canonical_product_name = art::canonicalProductName( + art::friendlyname::friendlyName("art::Assns"), + fMCTruthMCParticleAssnsInputModuleLabels.at(0).label(), + fMCTruthMCParticleAssnsInputModuleLabels.at(0).instance(), + fMCTruthMCParticleAssnsInputModuleLabels.at(0).process() + ); + auto product_id = art::ProductID(canonical_product_name); + std::cout << "Equivalent ProductID: " << product_id << std::endl; + } + InitializeCollectionMap(fMCTruthMCParticleAssnsInputModuleLabels,fMCTruthMCParticleAssnsMap); + std::cout << "fMCTruthMCParticleAssnsMap.size()=" << fMCTruthMCParticleAssnsMap.size() << std::endl; + + if(fVerbosity>1) + std::cout << "Finished initializing collection maps." << std::endl; +} + +template +void mix::SimInfoOverlayFilter::PutCollectionOntoEvent(art::Event & e, + CollectionMap & colmap) +{ + for(auto & cp : colmap) + e.put(std::move(cp.second),cp.first); +} + +void mix::SimInfoOverlayFilter::PutCollectionsOntoEvent(art::Event & e) +{ + PutCollectionOntoEvent(e,fMCTruthMap); + PutCollectionOntoEvent(e,fMCFluxMap); + PutCollectionOntoEvent(e,fGTruthMap); + PutCollectionOntoEvent(e,fBeamGateInfoMap); + + PutCollectionOntoEvent(e,fMCTruthMCFluxAssnsMap); + PutCollectionOntoEvent(e,fMCTruthGTruthAssnsMap); + + PutCollectionOntoEvent(e,fMCParticleMap); + PutCollectionOntoEvent(e,fSimEnergyDepositMap); + PutCollectionOntoEvent(e,fAuxDetSimChannelMap); + PutCollectionOntoEvent(e,fMCParticleLiteMap); + PutCollectionOntoEvent(e,fSimChannelMap); + PutCollectionOntoEvent(e,fSimPhotonsMap); + + PutCollectionOntoEvent(e,fMCTruthMCParticleAssnsMap); + + if(fVerbosity>1) + std::cout << "Finished putting collections onto the event." << std::endl; +} + +template +std::map< std::pair , art::Ptr > +mix::SimInfoOverlayFilter::FillCollectionMap(std::vector const& labels, + CollectionMap< std::vector > & colmap, + std::string type_name, + art::Event & e) +{ + std::map< std::pair , art::Ptr > + orig_artptr_lookup; + + for(auto label : labels){ + auto & out_ptrvec = colmap[label.instance()]; + art::PtrMaker makeArtPtr(e, label.instance()); + + gallery::Handle< std::vector > handle; + if(!gEvent.getByLabel< std::vector >(label,handle)) continue; + + //in the future, this will be the way we do it + //auto product_id = handle.id(); + + //but, that's only available in gallery v1_10_00 and up. + //For now, we gotta do it ourselves! + if(label.process().size()==0) + throw cet::exception("SimInfoOverlayFilter") + << "ERROR! All InputTags must be FULLY specified until gallery v1_10_00 is available." + << " InputTag: " << label << " not fully specified!"; + auto canonical_product_name = art::canonicalProductName(art::friendlyname::friendlyName(type_name), + label.label(), + label.instance(), + label.process()); + auto product_id = art::ProductID(canonical_product_name); + + if(fVerbosity>2) + std::cout << "Reading gallery handle. Canonical name = " << canonical_product_name << std::endl; + + auto const& vec(*handle); + for(size_t i_obj=0; i_objemplace_back(vec[i_obj]); + orig_artptr_lookup[std::make_pair(product_id,i_obj)] = makeArtPtr(out_ptrvec->size()-1); + } + } + + if(fVerbosity>2) + std::cout << "\tArtPtr lookup map has " << orig_artptr_lookup.size() << " entries." << std::endl; + + return orig_artptr_lookup; + +} + +template< class T, class U > +void mix::SimInfoOverlayFilter::FillAssnsCollectionMap(std::vector const& labels, + CollectionMap< art::Assns > & colmap, + std::map< std::pair , art::Ptr > const& + orig_artptr_lookup_left, + std::map< std::pair , art::Ptr > const& + orig_artptr_lookup_right) +{ + + for(auto label : labels){ + auto & out_ptrAssns = colmap[label.instance()]; + + gallery::Handle< art::Assns > handle; + if(!gEvent.getByLabel< art::Assns >(label,handle)) continue; + + auto const& assns(*handle); + for(size_t i_assn=0; i_assnaddSingle(newptr_left,newptr_right); + } + } + +} + +template< class T, class U, class D> +void mix::SimInfoOverlayFilter::FillAssnsCollectionMap(std::vector const& labels, + CollectionMap< art::Assns > & colmap, + std::map< std::pair , art::Ptr > const& + orig_artptr_lookup_left, + std::map< std::pair , art::Ptr > const& + orig_artptr_lookup_right) +{ + + std::cout << "Starting FillAssnsCollectionMap" << std::endl; + + // Notes on the inputs: + // colmap = fMCTruthMCParticleAssnsMap, CollectionMap filled with art::Assns + // orig_artptr_lookup_left, mctruth_artptr_lookup, map between product id and art::Ptr filled before calling this function + // orig_artptr_lookup_right, mcpart_artptr_lookup, map between product id and art::Ptr filled at the start of filter + + for(auto label : labels){ + std::cout << "label=" << label << std::endl; + auto & out_ptrAssns = colmap[label.instance()]; + + auto canonical_product_name = art::canonicalProductName( + art::friendlyname::friendlyName("art::Assns"), + label.label(), + label.instance(), + label.process() + ); + auto product_id = art::ProductID(canonical_product_name); + std::cout << "Equivalent ProductID: " << product_id << std::endl; + + gallery::Handle< art::Assns > handle; + if(!gEvent.getByLabel< art::Assns >(label,handle)) continue; + + auto const& assns(*handle); + for(size_t i_assn=0; i_assnaddSingle(newptr_left,newptr_right,data); + } + + } + +} + +bool mix::SimInfoOverlayFilter::filter(art::Event & e) +{ + + std::cout << std::endl << "STARTING FILTER" << std::endl; + + InitializeCollectionMaps(); + + //check if we have exhausted our simulation events. If so, we return false. + if(gEvent.atEnd()) { + if(fVerbosity>0) + std::cout << "We've reached the end of the simulation event stream. Returning false." << std::endl; + + PutCollectionsOntoEvent(e); + return false; + } + + if(fVerbosity>1){ + std::cout << "Processing input event: " + << "Run " << e.run() << ", " + << "Event " << e.event() << std::endl; + std::cout << "Processing simulation event: " + << "Run " << gEvent.eventAuxiliary().run() << ", " + << "Event " << gEvent.eventAuxiliary().event() << std::endl; + } + + if(fFillPOTInfo){ + /* + auto eventsInGalleryFile = gEvent.numberOfEventsInFile(); + gallery::Handle< sumdata::POTSummary > potsum_handle; + if(!gEvent.getByLabel(fPOTSummaryTag,potsum_handle)) + throw cet::exception("SimInfoOverlayFilter") << "No POTSummary object with tag " << fPOTSummaryTag; + + auto const& potsum(*potsum_handle); + fPOTSum_totpot += potsum.totpot/eventsInGalleryFile; + fPOTSum_totgoodpot += potsum.totgoodpot/eventsInGalleryFile; + fPOTSum_totspills += double(potsum.totspills)/eventsInGalleryFile; + fPOTSum_goodspills += double(potsum.goodspills)/eventsInGalleryFile; + */ + + fPOTSum_totpot += fSR_POTPerEvent[gEvent.eventAuxiliary().subRun()]; + fPOTSum_totgoodpot += fSR_GoodPOTPerEvent[gEvent.eventAuxiliary().subRun()]; + fPOTSum_totspills += fSR_SpillsPerEvent[gEvent.eventAuxiliary().subRun()]; + fPOTSum_goodspills += fSR_GoodSpillsPerEvent[gEvent.eventAuxiliary().subRun()]; + + } + + auto mctruth_artptr_lookup = FillCollectionMap(fMCTruthInputModuleLabels, + fMCTruthMap, + "std::vector", + e); + auto mcflux_artptr_lookup = FillCollectionMap(fMCFluxInputModuleLabels, + fMCFluxMap, + "std::vector", + e); + auto gtruth_artptr_lookup = FillCollectionMap(fGTruthInputModuleLabels, + fGTruthMap, + "std::vector", + e); + FillCollectionMap(fBeamGateInputModuleLabels, + fBeamGateInfoMap, + "std::vector", + e); + + FillAssnsCollectionMap(fMCTruthMCFluxAssnsInputModuleLabels, + fMCTruthMCFluxAssnsMap, + mctruth_artptr_lookup, + mcflux_artptr_lookup); + + FillAssnsCollectionMap(fMCTruthGTruthAssnsInputModuleLabels, + fMCTruthGTruthAssnsMap, + mctruth_artptr_lookup, + gtruth_artptr_lookup); + + std::cout << "Filling mcpart_artptr_lookup" << std::endl; + std::cout << "fMCParticleInputModuleLabels.size()= " << fMCParticleInputModuleLabels.size() << std::endl; + if(fMCParticleInputModuleLabels.size()){ + std::cout << "Using label " << fMCParticleInputModuleLabels.at(0) << std::endl; + auto canonical_product_name = art::canonicalProductName( + art::friendlyname::friendlyName("std::vector"), + fMCParticleInputModuleLabels.at(0).label(), + fMCParticleInputModuleLabels.at(0).instance(), + fMCParticleInputModuleLabels.at(0).process() + ); + auto product_id = art::ProductID(canonical_product_name); + std::cout << "Equivalent ProductID: " << product_id << std::endl; + } + + auto mcpart_artptr_lookup = FillCollectionMap(fMCParticleInputModuleLabels, + fMCParticleMap, + "std::vector", + e); + std::cout << "mcpart_artptr_lookup.size()=" << mcpart_artptr_lookup.size() << std::endl; + + + FillCollectionMap(fSimEnergyDepositInputModuleLabels, + fSimEnergyDepositMap, + "std::vector", + e); + + FillCollectionMap(fAuxDetSimChannelInputModuleLabels, + fAuxDetSimChannelMap, + "std::vector", + e); + + FillCollectionMap(fMCParticleLiteInputModuleLabels, + fMCParticleLiteMap, + "std::vector", + e); + + FillCollectionMap(fSimChannelInputModuleLabels, + fSimChannelMap, + "std::vector", + e); + + FillCollectionMap(fSimPhotonsInputModuleLabels, + fSimPhotonsMap, + "std::vector", + e); + + // this could be empty if running a "downstream" SimInfoMixer instance + if(fMCTruthInputModuleLabels.empty() && !fMCTruthMCParticleAssnsInputModuleLabels.empty()) { + // assume "generator" module generated the mctruth already + auto mclists = e.getMany>(); + for(size_t mcl = 0; mcl < mclists.size(); ++mcl){ + art::Handle< std::vector > mclistHandle = mclists[mcl]; + std::string process_name = mclistHandle.provenance()->processName(); + std::string instance_name = mclistHandle.provenance()->productInstanceName(); + std::string module_name = mclistHandle.provenance()->moduleLabel(); + // make sure we are only using products added in this session + if(process_name != moduleDescription().processName()) { + continue; + } + std::cout << "Making mctruth_artptr_lookup" << std::endl; + std::cout << "Using label " << module_name << " " << instance_name << " " << fMCTruthMCParticleAssnsInputModuleLabels.begin()->process() << std::endl; + auto canonical_product_name = art::canonicalProductName( + art::friendlyname::friendlyName("std::vector"), + module_name, + instance_name, + fMCTruthMCParticleAssnsInputModuleLabels.begin()->process()); + auto product_id = art::ProductID(canonical_product_name); + + // If the auxillary file was made with separate Gen/G4 steps - the wrong productID is assigned here + // causing a crash when tying to construct the assns, extra fhicl parameter to give the correct label + if(fAuxHasSeparateGenG4){ + std::cout << "Config indicates separate gen/g4 passes used in the aux file" << std::endl; + canonical_product_name = art::canonicalProductName( + art::friendlyname::friendlyName("std::vector"), + fMCTruthMCParticleAssnsMCTruthLookupLabel.at(0).label(), + fMCTruthMCParticleAssnsMCTruthLookupLabel.at(0).instance(), + fMCTruthMCParticleAssnsMCTruthLookupLabel.at(0).process()); + product_id = art::ProductID(canonical_product_name); + } + std::cout << "ProductID: " << product_id << std::endl; + for(size_t m = 0; m < mclistHandle->size(); ++m){ + art::Ptr mct(mclistHandle, m); + + mctruth_artptr_lookup[std::make_pair(product_id,m)] = mct; + } + } + } + + std::cout << "List of ProductIDs in mctruth_artptr_lookup:" << std::endl; + std::map< std::pair , art::Ptr >::iterator it_T; + for(it_T = mctruth_artptr_lookup.begin();it_T != mctruth_artptr_lookup.end();it_T++) + std::cout << it_T->first.first << std::endl; + + std::cout << "List of ProductIDs in mcpart_artptr_lookup:" << std::endl; + std::map< std::pair , art::Ptr >::iterator it_U; + for(it_U = mcpart_artptr_lookup.begin();it_U != mcpart_artptr_lookup.end();it_U++){ + std::cout << it_U->first.first << std::endl; + break; + } + + std::cout << "Calling FillAssnsCollectionMap" << std::endl; + FillAssnsCollectionMap + (fMCTruthMCParticleAssnsInputModuleLabels, + fMCTruthMCParticleAssnsMap, + mctruth_artptr_lookup, + mcpart_artptr_lookup); + std::cout << "Finished FillAssnsCollectionMap" << std::endl; + + //put onto event and loop the gallery event + PutCollectionsOntoEvent(e); + gEvent.next(); + return true; +} + +// POT map creation updated to handle multiple subruns in auxillary event file +void mix::SimInfoOverlayFilter::MakePOTMap(){ + + gEvent.first(); + + // Reset everything + fSR_POTPerEvent.clear(); + fSR_GoodPOTPerEvent.clear(); + fSR_SpillsPerEvent.clear(); + fSR_GoodSpillsPerEvent.clear(); + + unsigned int current_sr = 0; + unsigned int events_in_sr = 0; + + double current_sr_POT = 0; + double current_sr_GoodPOT = 0; + int current_sr_Spills = 0; + int current_sr_GoodSpills = 0; + + // Iterate through all events in auxillary file, count how many are in each subrun + while(!gEvent.atEnd()){ + + if(gEvent.eventAuxiliary().subRun() != current_sr){ + if(current_sr != 0){ + fSR_POTPerEvent[current_sr] = current_sr_POT/events_in_sr; + fSR_GoodPOTPerEvent[current_sr] = current_sr_GoodPOT/events_in_sr; + fSR_SpillsPerEvent[current_sr] = (double)current_sr_Spills/events_in_sr; + fSR_GoodSpillsPerEvent[current_sr] = (double)current_sr_GoodSpills/events_in_sr; + } + events_in_sr = 0; + current_sr = gEvent.eventAuxiliary().subRun(); + } + + gallery::Handle< sumdata::POTSummary > potsum_handle; + if(!gEvent.getByLabel(fPOTSummaryTag,potsum_handle)) + throw cet::exception("SimInfoOverlayFilterGenG4") << "No POTSummary object with tag " << fPOTSummaryTag; + + auto const& potsum(*potsum_handle); + current_sr_POT = potsum.totpot; + current_sr_GoodPOT = potsum.totgoodpot; + current_sr_Spills = potsum.totspills; + current_sr_GoodSpills = potsum.goodspills; + + events_in_sr++; + gEvent.next(); + } + + // Add the last sr + if(current_sr != 0){ + fSR_POTPerEvent[current_sr] = current_sr_POT/events_in_sr; + fSR_GoodPOTPerEvent[current_sr] = current_sr_GoodPOT/events_in_sr; + fSR_SpillsPerEvent[current_sr] = (double)current_sr_Spills/events_in_sr; + fSR_GoodSpillsPerEvent[current_sr] = (double)current_sr_GoodSpills/events_in_sr; + } + + // std::map::iterator it; + // for (it = fSR_POTPerEvent.begin(); it != fSR_POTPerEvent.end(); it++) std::cout << "SR=" << it->first << " POT/Event=" << it->second << std::endl; + + gEvent.first(); +} + +void mix::SimInfoOverlayFilter::beginJob() +{ +} + +bool mix::SimInfoOverlayFilter::beginRun(art::Run & r) +{ + return true; +} + +bool mix::SimInfoOverlayFilter::beginSubRun(art::SubRun & sr) +{ + if(fFillPOTInfo){ + fPOTSum_totpot = 0.0; + fPOTSum_totgoodpot = 0.0; + fPOTSum_totspills = 0.0; + fPOTSum_goodspills = 0.0; + } + return true; +} + +void mix::SimInfoOverlayFilter::endJob() +{ +} + +bool mix::SimInfoOverlayFilter::endRun(art::Run & r) +{ + return true; +} + +bool mix::SimInfoOverlayFilter::endSubRun(art::SubRun & sr) +{ + if(fFillPOTInfo){ + std::unique_ptr srpot_ptr(new sumdata::POTSummary()); + srpot_ptr->totpot = fPOTSum_totpot; + srpot_ptr->totgoodpot = fPOTSum_totgoodpot; + srpot_ptr->totspills = (int)fPOTSum_totspills; + srpot_ptr->goodspills = (int)fPOTSum_goodspills; + sr.put(std::move(srpot_ptr), art::subRunFragment()); + } + return true; +} + +void mix::SimInfoOverlayFilter::respondToCloseInputFile(art::FileBlock const & fb) +{ +} + +void mix::SimInfoOverlayFilter::respondToCloseOutputFiles(art::FileBlock const & fb) +{ +} + +void mix::SimInfoOverlayFilter::respondToOpenInputFile(art::FileBlock const &fb) +{ +} + +void mix::SimInfoOverlayFilter::respondToOpenOutputFiles(art::FileBlock const & fb) +{ +} + +DEFINE_ART_MODULE(mix::SimInfoOverlayFilter) \ No newline at end of file diff --git a/icaruscode/TPC/ICARUSWireCell/detsimmodules_wirecell_ICARUS.fcl b/icaruscode/TPC/ICARUSWireCell/detsimmodules_wirecell_ICARUS.fcl index fb6d9d486..d1430c9e3 100644 --- a/icaruscode/TPC/ICARUSWireCell/detsimmodules_wirecell_ICARUS.fcl +++ b/icaruscode/TPC/ICARUSWireCell/detsimmodules_wirecell_ICARUS.fcl @@ -54,6 +54,9 @@ icarus_simwire_wirecell: int_noise_scale: 1.0 coh_noise_scale: 1.0 + include_noise: true # noise on by default + overlay_drifter: false # default drifter by default + # Gain and shaping time //gain0: 14.9654 # mV/fC //gain1: 14.9654 # mV/fC diff --git a/icaruscode/TPC/ICARUSWireCell/icarus/sim.jsonnet b/icaruscode/TPC/ICARUSWireCell/icarus/sim.jsonnet index bd578ea57..9a7fb584c 100644 --- a/icaruscode/TPC/ICARUSWireCell/icarus/sim.jsonnet +++ b/icaruscode/TPC/ICARUSWireCell/icarus/sim.jsonnet @@ -88,5 +88,32 @@ function(params, tools) { signal: f.fanpipe('DepoSetFanout', self.signal_pipelines, 'FrameFanin', "simsignalgraph", outtags), splusn: f.fanpipe('DepoSetFanout', self.splusn_pipelines, 'FrameFanin', "simsplusngraph", outtags), + // Drifter for Overlay MC + overlay_drifter: g.pnode({ + local xregions = wc.unique_list(std.flattenArrays([v.faces for v in params.det.volumes])), + + type: "wclsICARUSDrifter", + data: params.lar { + rng: wc.tn(tools.random), + xregions: xregions, + time_offset: params.sim.depo_toffset, + + drift_speed: params.lar.drift_speed, + fluctuate: params.sim.fluctuate, + + DL: params.lar.DL, + DT: params.lar.DT, + lifetime: params.lar.lifetime, + ar39activity: 0, // no simulated activity + + // DB config + DBFileName: "tpc_elifetime_data", + DBTag: "v2r1", + ELifetimeCorrection: true, + Verbose: true, + TPC: 0, + }, + }, nin=1, nout=1, uses=[tools.random]), + } + sim, // tack on base for user sugar. }.ret diff --git a/icaruscode/TPC/ICARUSWireCell/icarus/wcls-multitpc-sim-drift-simchannel.jsonnet b/icaruscode/TPC/ICARUSWireCell/icarus/wcls-multitpc-sim-drift-simchannel.jsonnet index b5516b88f..9eb729ea3 100644 --- a/icaruscode/TPC/ICARUSWireCell/icarus/wcls-multitpc-sim-drift-simchannel.jsonnet +++ b/icaruscode/TPC/ICARUSWireCell/icarus/wcls-multitpc-sim-drift-simchannel.jsonnet @@ -168,11 +168,14 @@ local wcls_output = { }; //local deposio = io.numpy.depos(output); -local drifter = sim.drifter; + +local overlay_drifter = std.extVar("overlay_drifter"); + +local drifter = if overlay_drifter then sim.overlay_drifter else sim.drifter; local setdrifter = g.pnode({ type: 'DepoSetDrifter', data: { - drifter: "Drifter" + drifter: if overlay_drifter then "wclsICARUSDrifter" else "Drifter" } }, nin=1, nout=1, uses=[drifter]); @@ -331,6 +334,16 @@ local frame_summers = [ local actpipes = [g.pipeline([noises[n], coh_noises[n], digitizers[n], /*retaggers[n],*/ wcls_output.sim_digits[n]], name="noise-digitizer%d" %n) for n in std.range(0,3)]; local util = import 'pgrapher/experiment/icarus/funcs.jsonnet'; + +// local actpipes = [g.pipeline([noises[n], coh_noises[n], digitizers[n], /*retaggers[n],*/ wcls_output.sim_digits[n]], name="noise-digitizer%d" %n) for n in std.range(0,3)]; + +local include_noise = std.extVar("include_noise"); + +local actpipes = if include_noise then + [g.pipeline([noises[n], coh_noises[n], digitizers[n], wcls_output.h5io[n], wcls_output.sim_digits[n]], name="noise-digitizer%d" %n) for n in std.range(0,3)] + else + [g.pipeline([digitizers[n], wcls_output.h5io[n], wcls_output.sim_digits[n]], name="noise-digitizer%d" %n) for n in std.range(0,3)]; + local outtags = ['orig%d' % n for n in std.range(0, 3)]; local pipe_reducer = util.fansummer('DepoSetFanout', analog_pipes, frame_summers, actpipes, 'FrameFanin', 'fansummer', outtags); diff --git a/scripts/CMakeLists.txt b/scripts/CMakeLists.txt index c8c78a077..013cd33e7 100644 --- a/scripts/CMakeLists.txt +++ b/scripts/CMakeLists.txt @@ -1,2 +1,3 @@ install_scripts(SUBDIRS updates) +add_subdirectory(overlay_scripts) diff --git a/scripts/overlay_scripts/CMakeLists.txt b/scripts/overlay_scripts/CMakeLists.txt new file mode 100644 index 000000000..a8ce56a42 --- /dev/null +++ b/scripts/overlay_scripts/CMakeLists.txt @@ -0,0 +1 @@ +install_scripts() diff --git a/scripts/overlay_scripts/init_gen_common.sh b/scripts/overlay_scripts/init_gen_common.sh new file mode 100644 index 000000000..54e02dbaa --- /dev/null +++ b/scripts/overlay_scripts/init_gen_common.sh @@ -0,0 +1,30 @@ +#!/bin/sh + +MY_FCL_FILE=$1 #prodgenie_bnb_nu_filtered_NCPiZero_uboone.fcl +MY_N_EVENTS=$2 #2500 +MY_INPUT_DATA_FILE=$3 #data offbeamminbias file used for overlays +MY_OUTPUT_FILE=genfile.root.local +MY_OUT_LOG=larInitGen.out +MY_ERR_LOG=larInitGen.err + +echo "#include \"$MY_FCL_FILE\"" > local_gen.fcl +#echo "physics.producers.generator.FluxCopyMethod: \"IFDH\"" >> local_gen.fcl +#echo "physics.producers.generator.MaxFluxFileMB: 500" >> local_gen.fcl +echo "services.IFDH: {}" >> local_gen.fcl + +echo "#include \"$MY_FCL_FILE\"" > local_gen_include.fcl +#echo "physics.producers.generator.FluxCopyMethod: \"IFDH\"" >> local_gen_include.fcl +#echo "physics.producers.generator.MaxFluxFileMB: 500" >> local_gen_include.fcl +echo "services.IFDH: {}" >> local_gen_include.fcl +echo "gen_detail: { physics: {@table::physics} services: {@table::services} outputs: {@table::outputs} source: {@table::source} process_name: @local::process_name }" >> local_gen_include.fcl + +if [ -f $MY_OUTPUT_FILE ]; then + echo "File $MY_OUTPUT_FILE exists, so not generating again." +else + echo "File $MY_OUTOUT_FILE does not exist, so running generation of it..." + echo "Running command 'lar -c local_gen.fcl -T ./genfile_hist.root -o $MY_OUTPUT_FILE -n $MY_N_EVENTS > $MY_OUT_LOG 2> $MY_ERR_LOG' " + ##lar -c local_gen.fcl -T ./genfile_hist.root -o genfile.root.temp -n $MY_N_EVENTS + echo "Running command 'lar -c intime_gen_overlay_SimInfoMixer1.fcl -s $MY_INPUT_DATA_FILE -T ./genfile_pot_hist.root -n -1' " + echo "Running command 'lar -c intime_gen_overlay_SimInfoMixer2.fcl -s *_simmxd1.root -T ./genfile_pot_hist.root -n -1' " + echo "Running command 'rm genfile.root.temp' " +fi diff --git a/scripts/overlay_scripts/init_gen_intime.sh b/scripts/overlay_scripts/init_gen_intime.sh new file mode 100644 index 000000000..364c81604 --- /dev/null +++ b/scripts/overlay_scripts/init_gen_intime.sh @@ -0,0 +1,94 @@ +#!/bin/sh + +############################################################ +# Help # +############################################################ +Help() +{ + # Display Help + echo "To run keepup_def_making.sh" + echo + echo "Syntax: sh keepup_def_making.sh [-s -v -f -l -m | -h ]" + echo "options:" + echo "s Enter stage name. Default is 'raw'." + echo "v Enter software version. Default is 'v09_56_00_01'" + echo "f Enter starting run number for range. Default is '9301'" + echo "l Enter ending run number for range. Deafult is '9500'" + echo "m Enter 'y' to make definitions. Default is 'n'" + echo "h Print this Help menu." + echo +} + +############################################################ +# Process the input options. # +############################################################ + + +FCLFILE="prodcorsika_proton_intime_icarus_bnb_sce_on.fcl" +NEVT=50 +GNEVT=2500 +GENFILE="genfile.root.local" +OUTPUTFILE="gen_1234.root" +OUTGENFILE="mcgen.root" +INPUTDATAFILE="/path/to/file.root" + +MY_OUT_LOG=larInitGen.out +MY_ERR_LOG=larInitGen.err + +# Get the options +while getopts ":h:s:n:g:o:s:c:f:" option; do + case $option in + h) # display Help + Help + exit;; + s) # Enter input data file + INPUTDATAFILE=$OPTARG;; + n) # Enter number of events to process + NEVT=$OPTARG;; + g) #Enter number of events to generate + GNEVT=$OPTARG;; + o) #Enter the output file name + OUTPUTFILE=$OPTARG;; + s) #name of input data file + INPUTDATAFILE=$OPTARG;; + c) #Name of input fcl file + FCLFILE=$OPTARG;; + f) #Name of output genfile + OUTGENFILE=$OPTARG;; + \?) # Invalid option + echo "Error: Invalid option" + exit;; + esac +done + + + +echo "#include \"$FCLFILE\"" > local_gen.fcl +#echo "physics.producers.generator.FluxCopyMethod: \"IFDH\"" >> local_gen.fcl +#echo "physics.producers.generator.MaxFluxFileMB: 500" >> local_gen.fcl +echo "services.IFDH: {}" >> local_gen.fcl +cat local_gen.fcl + +echo "#include \"$FCLFILE\"" > local_gen_include.fcl +#echo "physics.producers.generator.FluxCopyMethod: \"IFDH\"" >> local_gen_include.fcl +#echo "physics.producers.generator.MaxFluxFileMB: 500" >> local_gen_include.fcl +echo "services.IFDH: {}" >> local_gen_include.fcl +echo "gen_detail: { physics: {@table::physics} services: {@table::services} outputs: {@table::outputs} source: {@table::source} process_name: @local::process_name }" >> local_gen_include.fcl + + + +if [ -f $GENFILE ]; then + echo "File $GENFILE exists, so not generating again." +else + echo "File $GENFILE does not exist, so running generation of it..." + echo "Running command 'lar -c local_gen.fcl -T ./genfile_hist.root -o $GENFILE -n $GNEVT > $MY_OUT_LOG 2> $MY_ERR_LOG' " + lar -c local_gen.fcl -T ./genfile_hist.root -o $GENFILE -n $GNEVT > $MY_OUT_LOG 2> $MY_ERR_LOG +fi +echo "fcl_file: $FCLFILE; GNEVT: $GNEVT; NEVT: $NEVT; OUTPUTFILE: $OUTPUTFILE; INPUTDATAFILE: $INPUTDATAFILE" + +echo "Running command 'lar -c intime_gen_overlay_SimInfoMixer1.fcl -s $INPUTDATAFILE -T ./genfile_pot_hist.root -n -1' " +echo "Running command 'lar -c intime_gen_overlay_SimInfoMixer2.fcl -s *_simmxd1.root -o $OUTPUTFILE -T ./genfile_pot_hist.root -n -1' " +lar -c intime_gen_overlay_SimInfoMixer1.fcl -s $INPUTDATAFILE -T ./genfile_pot_hist.root -n $NEVT +lar -c intime_gen_overlay_SimInfoMixer2.fcl -s *_simmxd1.root -o $OUTPUTFILE -T ./genfile_pot_hist.root -n $NEVT +cp $GENFILE $OUTGENFILE +rm *_simmxd1.root