LArSoft
Software for Liquid Argon time projection chambers
Explicit code changes since v06_18_00
These are changes in larsoft v06_18_01.
larwirecell
larwirecell/LArInterface/Test_module.cc
diff --git a/larwirecell/LArInterface/Test_module.cc b/larwirecell/LArInterface/Test_module.cc
deleted file mode 100644
index ce609e9..0000000
--- a/larwirecell/LArInterface/Test_module.cc
+++ /dev/null
@@ -1,243 +0,0 @@
-#ifndef TESTMODULE_H
-#define TESTMODULE_H
-
-#include <string>
-#include <vector>
-#include <iostream>
-
-// Framework includes
-#include "art/Framework/Principal/Event.h"
-//#include "fhiclcpp/ParameterSet.h"
-
-#include "art/Framework/Core/ModuleMacros.h"
-#include "art/Framework/Core/EDProducer.h"
-#include "art/Framework/Principal/Handle.h"
-#include "art/Framework/Services/Optional/TFileService.h"
-#include "art/Framework/Services/Registry/ServiceHandle.h"
-#include "messagefacility/MessageLogger/MessageLogger.h"
-
-#include "lardataobj/RawData/RawDigit.h"
-#include "lardataobj/RawData/raw.h"
-
-#include <numeric> // iota
-#include "WireCellIface/SimpleFrame.h"
-#include "WireCellIface/SimpleTrace.h"
-#include "WireCellSigProc/OmnibusNoiseFilter.h"
-#include "WireCellSigProc/OneChannelNoise.h"
-#include "WireCellSigProc/CoherentNoiseSub.h"
-#include "WireCellSigProc/SimpleChannelNoiseDB.h"
-
-using namespace WireCell;
-using namespace std;
-
-//class TestTrace;
-
-namespace testalg {
-
- class Test : public art::EDProducer {
-
- public:
- explicit Test(fhicl::ParameterSet const& pset);
- virtual ~Test();
-
- void produce(art::Event & evt);
- void reconfigure(fhicl::ParameterSet const& pset);
-
- void beginJob();
- void endJob();
-
- private:
-
- //******************************
- //Variables Taken from FHICL File
- //std::string fRawDigitModuleLabel; //label for rawdigit module
- }; //end class Noise
-
-
- //-------------------------------------------------------------------
- Test::Test(fhicl::ParameterSet const& pset)
- : EDProducer(){
- this->reconfigure(pset);
- produces<std::vector<raw::RawDigit> >();
- }
-
- //-------------------------------------------------------------------
- Test::~Test(){}
-
- //-------------------------------------------------------------------
- void Test::reconfigure(fhicl::ParameterSet const& pset){
- }
-
- //-------------------------------------------------------------------
- void Test::beginJob(){
- art::ServiceHandle<art::TFileService> tfs;
- //art::ServiceHandle<util::LArWireCellNoiseFilterService> larWireCellNF;
- //larWireCellNF->print(17);
- }
-
- //-------------------------------------------------------------------
- void Test::endJob(){
- art::ServiceHandle<art::TFileService> tfs;
- }
-
- //-------------------------------------------------------------------
- void Test::produce(art::Event & evt){
- art::Handle< std::vector<raw::RawDigit> > rawDigitHandle;
- evt.getByLabel("daq",rawDigitHandle);
- std::vector<raw::RawDigit> const& rawDigitVector(*rawDigitHandle);
- const unsigned int n_channels = rawDigitVector.size();
-
- //skip empty events
- if( n_channels == 0 )
- return;
-
- // S&C microboone sampling parameter database
- const double tick = 0.5*units::microsecond;
- const int nsamples = 9595;
-
- // Q&D microboone channel map
- std::vector<int> uchans(2400), vchans(2400), wchans(3456);
- const int nchans = uchans.size() + vchans.size() + wchans.size();
- std::iota(uchans.begin(), uchans.end(), 0);
- std::iota(vchans.begin(), vchans.end(), vchans.size());
- std::iota(wchans.begin(), wchans.end(), vchans.size() + uchans.size());
-
- // Q&D nominal baseline
- const double unombl=2048.0, vnombl=2048.0, wnombl=400.0;
-
- // Q&D miss-configured channel database
- vector<int> miscfgchan;
- const double from_gain_mVfC=7.8, to_gain_mVfC=14.0,
- from_shaping=1.0*units::microsecond, to_shaping=2.0*units::microsecond;
- for (int ind=2016; ind<= 2095; ++ind) { miscfgchan.push_back(ind); }
- for (int ind=2192; ind<= 2303; ++ind) { miscfgchan.push_back(ind); }
- for (int ind=2352; ind< 2400; ++ind) { miscfgchan.push_back(ind); }
-
- // hard-coded bad channels
- vector<int> bad_channels;
- for (unsigned int i=0;i!=wchans.size();i++){
- if (i>=7136 - 4800 && i <=7263 - 4800){
- if (i != 7200- 4800 && i!=7215 - 4800)
- bad_channels.push_back(i+4800);
- }
- }
-
- // Q&D RC+RC time constant - all have same.
- const double rcrc = 1.0*units::millisecond;
- vector<int> rcrcchans(nchans);
- std::iota(rcrcchans.begin(), rcrcchans.end(), 0);
-
- //harmonic noises
- vector<int> harmonicchans(uchans.size() + vchans.size());
- std::iota(harmonicchans.begin(), harmonicchans.end(), 0);
-
- vector<int> special_chans;
- special_chans.push_back(2240);
-
- WireCellSigProc::SimpleChannelNoiseDB::mask_t h36kHz(0,169,173);
- WireCellSigProc::SimpleChannelNoiseDB::mask_t h108kHz(0,513,516);
- WireCellSigProc::SimpleChannelNoiseDB::mask_t hspkHz(0,17,19);
- WireCellSigProc::SimpleChannelNoiseDB::multimask_t hharmonic;
- hharmonic.push_back(h36kHz);
- hharmonic.push_back(h108kHz);
- WireCellSigProc::SimpleChannelNoiseDB::multimask_t hspecial;
- hspecial.push_back(h36kHz);
- hspecial.push_back(h108kHz);
- hspecial.push_back(hspkHz);
-
- // do the coherent subtraction
- std::vector< std::vector<int> > channel_groups;
- for (unsigned int i=0;i!=172;i++){
- //for (int i=150;i!=151;i++){
- std::vector<int> channel_group;
- for (int j=0;j!=48;j++){
- channel_group.push_back(i*48+j);
- }
- channel_groups.push_back(channel_group);
- }
-
- auto noise = new WireCellSigProc::SimpleChannelNoiseDB;
- // initialize
- noise->set_sampling(tick, nsamples);
- // set nominal baseline
- noise->set_nominal_baseline(uchans, unombl);
- noise->set_nominal_baseline(vchans, vnombl);
- noise->set_nominal_baseline(wchans, wnombl);
- // set misconfigured channels
- noise->set_gains_shapings(miscfgchan, from_gain_mVfC, to_gain_mVfC, from_shaping, to_shaping);
- // do the RCRC
- noise->set_rcrc_constant(rcrcchans, rcrc);
- // set initial bad channels
- noise->set_bad_channels(bad_channels);
- // set the harmonic filter
- noise->set_filter(harmonicchans,hharmonic);
- noise->set_filter(special_chans,hspecial);
- noise->set_channel_groups(channel_groups);
- //Define database object
- shared_ptr<WireCell::IChannelNoiseDatabase> noise_sp(noise);
-
- auto one = new WireCellSigProc::OneChannelNoise;
- one->set_channel_noisedb(noise_sp);
- shared_ptr<WireCell::IChannelFilter> one_sp(one);
- auto many = new WireCellSigProc::CoherentNoiseSub;
- shared_ptr<WireCell::IChannelFilter> many_sp(many);
-
- //define noisefilter object
- WireCellSigProc::OmnibusNoiseFilter bus;
- bus.set_channel_filters({one_sp});
- bus.set_grouped_filters({many_sp});
- bus.set_channel_noisedb(noise_sp);
-
- //load waveforms into traces
- ITrace::vector traces;
- for(unsigned int ich=0; ich<n_channels; ich++){
- const size_t n_samp = rawDigitVector.at(ich).NADC();
- if( n_samp == 0 )
- continue;
-
- ITrace::ChargeSequence charges;
- for( unsigned int s = 0 ; s < n_samp ; s++ ){
- float q = (float)rawDigitVector.at(ich).ADCs().at(s);
- charges.push_back(q);
- }
- unsigned int chan = rawDigitVector.at(ich).Channel();
- WireCell::SimpleTrace* st = new WireCell::SimpleTrace(chan, 0.0, charges);
- traces.push_back(ITrace::pointer(st));
- }
-
- //Load traces into frame
- WireCell::SimpleFrame* sf = new WireCell::SimpleFrame(0, 0, traces);
- IFrame::pointer frame = IFrame::pointer(sf);
- IFrame::pointer quiet;
-
- //Do filtering
- bus(frame, quiet);
-
- //Output results
- std::unique_ptr<std::vector<raw::RawDigit> > filteredRawDigit(new std::vector<raw::RawDigit>);
- std::vector< short > waveform;
-
- auto quiet_traces = quiet->traces();
- for (auto quiet_trace : *quiet_traces.get()) {
- //int tbin = quiet_trace->tbin();
- unsigned int ch = quiet_trace->channel();
- auto quiet_charges = quiet_trace->charge();
- int counter = 0;
- waveform.clear();
- for (auto q : quiet_charges) {
- waveform.push_back(q);
- counter++;
- }
- unsigned int n_samp = waveform.size();
- filteredRawDigit->emplace_back( raw::RawDigit( ch , n_samp, waveform, raw::kNone) );
- }
-
- //filtered raw digits
- evt.put(std::move(filteredRawDigit));
- }
-
- DEFINE_ART_MODULE(Test)
-
-} //end namespace testalg
-
-#endif //TESTMODULE_H
larwirecell/LArInterface/WireCellNoiseFilter_module.cc
diff --git a/larwirecell/LArInterface/WireCellNoiseFilter_module.cc b/larwirecell/LArInterface/WireCellNoiseFilter_module.cc
new file mode 100644
index 0000000..ae36738
--- /dev/null
+++ b/larwirecell/LArInterface/WireCellNoiseFilter_module.cc
@@ -0,0 +1,323 @@
+#ifndef WIRECELLNOISEFILTERMODULE_H // <--- this is not necessary since this file is not included...
+#define WIRECELLNOISEFILTERMODULE_H
+
+// Framework includes
+#include "art/Framework/Principal/Event.h"
+
+#include "art/Framework/Core/ModuleMacros.h"
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Services/Optional/TFileService.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "lardataobj/RawData/RawDigit.h"
+#include "lardataobj/RawData/raw.h"
+
+#include <numeric> // iota
+#include "WireCellUtil/Units.h"
+#include "WireCellIface/SimpleFrame.h"
+#include "WireCellIface/SimpleTrace.h"
+#include "WireCellSigProc/OmnibusNoiseFilter.h"
+#include "WireCellSigProc/OneChannelNoise.h"
+#include "WireCellSigProc/CoherentNoiseSub.h"
+#include "WireCellSigProc/SimpleChannelNoiseDB.h"
+
+// Undefine the def in WireCellUtil/Units.h so we can regain the CLHEP definitions in the below header files
+// They don't have the "units" namespace...
+#undef HEP_SYSTEM_OF_UNITS_H
+
+#include "larevt/CalibrationDBI/Interface/DetPedestalService.h"
+#include "larevt/CalibrationDBI/Interface/DetPedestalProvider.h"
+#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
+#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
+#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
+#include "larcore/Geometry/Geometry.h"
+
+#include <string>
+#include <vector>
+#include <iostream>
+
+using namespace WireCell;
+//using namespace std;
+
+namespace noisefilteralg {
+
+class WireCellNoiseFilter : public art::EDProducer {
+
+public:
+ explicit WireCellNoiseFilter(fhicl::ParameterSet const& pset);
+ virtual ~WireCellNoiseFilter();
+
+ void produce(art::Event & evt);
+ void reconfigure(fhicl::ParameterSet const& pset);
+
+ void beginJob();
+ void endJob();
+
+private:
+
+ void DoNoiseFilter(const std::vector<raw::RawDigit>&, std::vector<raw::RawDigit>&) const;
+
+ //******************************
+ //Variables Taken from FHICL File
+ std::string fDigitModuleLabel; //label for rawdigit module
+ bool fDoNoiseFiltering;
+ bool fTruncateTicks;
+ size_t fWindowSize;
+ size_t fNumTicksToDropFront;
+
+ // services
+}; //end class Noise
+
+
+//-------------------------------------------------------------------
+WireCellNoiseFilter::WireCellNoiseFilter(fhicl::ParameterSet const& pset)
+ : EDProducer()
+{
+ this->reconfigure(pset);
+ produces<std::vector<raw::RawDigit> >();
+}
+
+//-------------------------------------------------------------------
+WireCellNoiseFilter::~WireCellNoiseFilter(){}
+
+//-------------------------------------------------------------------
+void WireCellNoiseFilter::reconfigure(fhicl::ParameterSet const& pset){
+ fDigitModuleLabel = pset.get<std::string>("DigitModuleLabel", "daq");
+ fDoNoiseFiltering = pset.get<bool> ("DoNoiseFiltering", true );
+ fTruncateTicks = pset.get<bool> ("TruncateTicks", true );
+ fWindowSize = pset.get<size_t> ("WindowSize", 6400 );
+ fNumTicksToDropFront = pset.get<size_t> ("NumTicksToDropFront", 2400 );
+}
+
+//-------------------------------------------------------------------
+void WireCellNoiseFilter::beginJob(){
+ art::ServiceHandle<art::TFileService> tfs;
+ //art::ServiceHandle<util::LArWireCellNoiseFilterService> larWireCellNF;
+ //larWireCellNF->print(17);
+}
+
+//-------------------------------------------------------------------
+void WireCellNoiseFilter::endJob(){
+ art::ServiceHandle<art::TFileService> tfs;
+}
+
+//-------------------------------------------------------------------
+void WireCellNoiseFilter::produce(art::Event & evt)
+{
+ // Recover services we will need
+ const lariov::DetPedestalProvider& pedestalValues = art::ServiceHandle<lariov::DetPedestalService>()->GetPedestalProvider();
+ const detinfo::DetectorProperties& detectorProperties = *lar::providerFrom<detinfo::DetectorPropertiesService>();
+ larana/T0Finder/PhotonCounterT0Matching_module.cc
+ art::Handle< std::vector<raw::RawDigit> > rawDigitHandle;
+ evt.getByLabel(fDigitModuleLabel,rawDigitHandle);
+
+ // Define the output vector (in case we don't do anything)
+ std::unique_ptr<std::vector<raw::RawDigit> > filteredRawDigit(new std::vector<raw::RawDigit>);
+
+ if (rawDigitHandle.isValid() && rawDigitHandle->size() > 0)
+ {
+ const std::vector<raw::RawDigit>& rawDigitVector(*rawDigitHandle);
+
+ if (fDoNoiseFiltering) DoNoiseFilter(rawDigitVector, *filteredRawDigit);
+ else
+ {
+ // Enable truncation
+ size_t startBin(0);
+ size_t stopBin(detectorProperties.NumberTimeSamples());
+
+ if (fTruncateTicks)
+ {
+ startBin = fNumTicksToDropFront;
+ stopBin = fNumTicksToDropFront + fWindowSize;
+ }
+
+ raw::RawDigit::ADCvector_t outputVector(detectorProperties.NumberTimeSamples());
+
+ for(const auto& rawDigit : rawDigitVector)
+ {
+ if (rawDigit.NADC() < detectorProperties.NumberTimeSamples()) continue;
+
+ const raw::RawDigit::ADCvector_t& rawAdcVec = rawDigit.ADCs();
+
+ unsigned int channel = rawDigit.Channel();
+ float pedestal = pedestalValues.PedMean(channel);
+
+ std::copy(rawAdcVec.begin() + startBin, rawAdcVec.begin() + stopBin, outputVector.begin());
+
+ filteredRawDigit->emplace_back( raw::RawDigit( channel , outputVector.size(), outputVector, raw::kNone) );
+ filteredRawDigit->back().SetPedestal(pedestal,2.0);
+ }
+ }
+ }
+
+ //filtered raw digits
+ evt.put(std::move(filteredRawDigit));
+}
+
+void WireCellNoiseFilter::DoNoiseFilter(const std::vector<raw::RawDigit>& inputWaveforms, std::vector<raw::RawDigit>& outputWaveforms) const
+{
+ // Recover services we will need
+ const lariov::ChannelStatusProvider& channelStatus = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
+ const lariov::DetPedestalProvider& pedestalValues = art::ServiceHandle<lariov::DetPedestalService>()->GetPedestalProvider();
+ const geo::GeometryCore& geometry = *lar::providerFrom<geo::Geometry>();
+ const detinfo::DetectorProperties& detectorProperties = *lar::providerFrom<detinfo::DetectorPropertiesService>();
+
+ const unsigned int n_channels = inputWaveforms.size();
+
+ // S&C microboone sampling parameter database
+ const double tick = detectorProperties.SamplingRate(); // 0.5 * units::microsecond;
+ const size_t nsamples = detectorProperties.NumberTimeSamples();
+
+ // Q&D microboone channel map
+ std::vector<int> uchans(geometry.Nwires(0)), vchans(geometry.Nwires(1)), wchans(geometry.Nwires(2));
+ const int nchans = uchans.size() + vchans.size() + wchans.size();
+ std::iota(uchans.begin(), uchans.end(), 0);
+ std::iota(vchans.begin(), vchans.end(), vchans.size());
+ std::iota(wchans.begin(), wchans.end(), vchans.size() + uchans.size());
+
+ // Q&D nominal baseline
+ const double unombl=2048.0, vnombl=2048.0, wnombl=400.0;
+
+ // Q&D miss-configured channel database
+ std::vector<int> miscfgchan;
+ const double from_gain_mVfC=7.8, to_gain_mVfC=14.0,from_shaping=1.0*units::microsecond, to_shaping=2.0*units::microsecond;
+ for (int ind=2016; ind<= 2095; ++ind) { miscfgchan.push_back(ind); }
+ for (int ind=2192; ind<= 2303; ++ind) { miscfgchan.push_back(ind); }
+ for (int ind=2352; ind< 2400; ++ind) { miscfgchan.push_back(ind); }
+
+ // Recover bad channels from the database
+ std::vector<int> bad_channels;
+ for(int channelIdx=0; channelIdx<nchans; channelIdx++) if (channelStatus.IsBad(channelIdx)) bad_channels.push_back(channelIdx);
+
+ // Q&D RC+RC time constant - all have same.
+ const double rcrc = 1.0*units::millisecond;
+ std::vector<int> rcrcchans(nchans);
+ std::iota(rcrcchans.begin(), rcrcchans.end(), 0);
+
+ //harmonic noises
+ std::vector<int> harmonicchans(uchans.size() + vchans.size());
+ std::iota(harmonicchans.begin(), harmonicchans.end(), 0);
+
+ std::vector<int> special_chans;
+ special_chans.push_back(2240);
+
+ WireCellSigProc::SimpleChannelNoiseDB::mask_t h36kHz(0,169,173);
+ WireCellSigProc::SimpleChannelNoiseDB::mask_t h108kHz(0,513,516);
+ WireCellSigProc::SimpleChannelNoiseDB::mask_t hspkHz(0,17,19);
+ WireCellSigProc::SimpleChannelNoiseDB::multimask_t hharmonic;
+ hharmonic.push_back(h36kHz);
+ hharmonic.push_back(h108kHz);
+ WireCellSigProc::SimpleChannelNoiseDB::multimask_t hspecial;
+ hspecial.push_back(h36kHz);
+ hspecial.push_back(h108kHz);
+ hspecial.push_back(hspkHz);
+
+ // do the coherent subtraction
+ std::vector< std::vector<int> > channel_groups;
+ for (unsigned int i=0;i!=172;i++){
+ //for (int i=150;i!=151;i++){
+ std::vector<int> channel_group;
+ for (int j=0;j!=48;j++){
+ channel_group.push_back(i*48+j);
+ }
+ channel_groups.push_back(channel_group);
+ }
+
+ auto noise = new WireCellSigProc::SimpleChannelNoiseDB;
+ // initialize
+ noise->set_sampling(tick, fWindowSize);
+ // set nominal baseline
+ noise->set_nominal_baseline(uchans, unombl);
+ noise->set_nominal_baseline(vchans, vnombl);
+ noise->set_nominal_baseline(wchans, wnombl);
+ // set misconfigured channels
+ noise->set_gains_shapings(miscfgchan, from_gain_mVfC, to_gain_mVfC, from_shaping, to_shaping);
+ // do the RCRC
+ noise->set_rcrc_constant(rcrcchans, rcrc);
+ // set initial bad channels
+ noise->set_bad_channels(bad_channels);
+ // set the harmonic filter
+ noise->set_filter(harmonicchans,hharmonic);
+ noise->set_filter(special_chans,hspecial);
+ noise->set_channel_groups(channel_groups);
+ //Define database object
+ std::shared_ptr<WireCell::IChannelNoiseDatabase> noise_sp(noise);
+
+ auto one = new WireCellSigProc::OneChannelNoise;
+ one->set_channel_noisedb(noise_sp);
+ std::shared_ptr<WireCell::IChannelFilter> one_sp(one);
+ auto many = new WireCellSigProc::CoherentNoiseSub;
+ std::shared_ptr<WireCell::IChannelFilter> many_sp(many);
+
+ //define noisefilter object
+ WireCellSigProc::OmnibusNoiseFilter bus;
+ bus.set_channel_filters({one_sp});
+ bus.set_grouped_filters({many_sp});
+ bus.set_channel_noisedb(noise_sp);
+
+ // Enable truncation
+ size_t startBin(0);
+ size_t stopBin(nsamples);
+
+ if (fTruncateTicks)
+ {
+ startBin = fNumTicksToDropFront;
+ stopBin = fNumTicksToDropFront + fWindowSize;
+ }
+
+ //load waveforms into traces
+ ITrace::vector traces;
+ for(unsigned int ich=0; ich<n_channels; ich++)
+ {
+ if( inputWaveforms.at(ich).NADC() < nsamples ) continue;
+
+ const raw::RawDigit::ADCvector_t& rawAdcVec = inputWaveforms.at(ich).ADCs();
+
+ ITrace::ChargeSequence charges;
+
+ charges.resize(fWindowSize);
+
+ std::transform(rawAdcVec.begin() + startBin, rawAdcVec.begin() + stopBin, charges.begin(), [](auto& adcVal){return float(adcVal);});
+
+ unsigned int chan = inputWaveforms.at(ich).Channel();
+ WireCell::SimpleTrace* st = new WireCell::SimpleTrace(chan, 0.0, charges);
+ traces.push_back(ITrace::pointer(st));
+ }
+
+ //Load traces into frame
+ WireCell::SimpleFrame* sf = new WireCell::SimpleFrame(0, 0, traces);
+ IFrame::pointer frame = IFrame::pointer(sf);
+ IFrame::pointer quiet;
+
+ //Do filtering
+ bus(frame, quiet);
+
+ //Output results
+ std::vector< short > waveform(fWindowSize);
+
+ auto quiet_traces = quiet->traces();
+ for (auto quiet_trace : *quiet_traces.get()) {
+ //int tbin = quiet_trace->tbin();
+ unsigned int channel = quiet_trace->channel();
+ auto quiet_charges = quiet_trace->charge();
+
+ // Recover the database version of the pedestal, we'll offset the waveforms so it matches
+ float pedestal = pedestalValues.PedMean(channel);
+
+ std::transform(quiet_charges.begin(), quiet_charges.end(), waveform.begin(), [pedestal](auto charge){return std::round(charge+pedestal);});
+
+ outputWaveforms.emplace_back( raw::RawDigit( channel , waveform.size(), waveform, raw::kNone) );
+ outputWaveforms.back().SetPedestal(pedestal,2.0);
+ }
+
+ return;
+}
+
+
+DEFINE_ART_MODULE(WireCellNoiseFilter)
+
+} //end namespace WireCellNoiseFilteralg
+
+#endif //WireCellNoiseFilterMODULE_H
larwirecell/LArInterface/wirecellmodules.fcl
diff --git a/larwirecell/LArInterface/wirecellmodules.fcl b/larwirecell/LArInterface/wirecellmodules.fcl
new file mode 100644
index 0000000..90af2f9
--- /dev/null
+++ b/larwirecell/LArInterface/wirecellmodules.fcl
@@ -0,0 +1,16 @@
+// from larreco/RecoAlg:
+#include "hitalgorithms.fcl"
+
+BEGIN_PROLOG
+
+wirecellnoisefilter:
+{
+ module_type: "WireCellNoiseFilter"
+ DigitModuleLabel: "daq" # RawDigit producer
+ DoNoiseFiltering: true # Filter noise, false = pass through
+ TruncateTicks: false # Truncate the waveoforms?
+ WindowSize: 9600 # Total size of waveform
+ NumTicksToDropFront: 0 # Number ticks to drop from front
+}
+
+END_PROLOG
larwirecell/LArInterface/wirecellmodules_microboone.fcl
diff --git a/larwirecell/LArInterface/wirecellmodules_microboone.fcl b/larwirecell/LArInterface/wirecellmodules_microboone.fcl
new file mode 100644
index 0000000..fad2929
--- /dev/null
+++ b/larwirecell/LArInterface/wirecellmodules_microboone.fcl
@@ -0,0 +1,11 @@
+#include "wirecellmodules.fcl"
+
+BEGIN_PROLOG
+
+microboone_wirecellnoisefilter: @local::wirecellnoisefilter
+microboone_wirecellnoisefilter.DoNoiseFiltering: true
+microboone_wirecellnoisefilter.TruncateTicks: true
+microboone_wirecellnoisefilter.WindowSize: 6400
+microboone_wirecellnoisefilter.NumTicksToDropFront: 2400
+
+END_PROLOG
larsim
larsim/EventGenerator/MARLEY/MARLEYGen_module.cc
diff --git a/larsim/EventGenerator/MARLEY/MARLEYGen_module.cc b/larsim/EventGenerator/MARLEY/MARLEYGen_module.cc
index 259f84b..6648eca 100644
--- a/larsim/EventGenerator/MARLEY/MARLEYGen_module.cc
+++ b/larsim/EventGenerator/MARLEY/MARLEYGen_module.cc
@@ -493,8 +493,8 @@ void evgen::MarleyGen::reconfigure_marley(const fhicl::ParameterSet& p)
fhicl::ParameterSet source_params = p.get<fhicl::ParameterSet>("source");
- prepare_reactions(reactions);
prepare_structure(structure);
+ prepare_reactions(reactions);
prepare_neutrino_source(source_params);
}
larsim/MCCheater/PhotonBackTracker.h
diff --git a/larsim/MCCheater/PhotonBackTracker.h b/larsim/MCCheater/PhotonBackTracker.h
index 28fe9ae..c222539 100644
--- a/larsim/MCCheater/PhotonBackTracker.h
+++ b/larsim/MCCheater/PhotonBackTracker.h
@@ -42,7 +42,7 @@ namespace cheat{
public:
PhotonBackTracker(fhicl::ParameterSet const& pset,
- art::ActivityRegistry& reg);
+ art::ActivityRegistry& reg);
~PhotonBackTracker();
void reconfigure(fhicl::ParameterSet const& pset);
@@ -80,7 +80,7 @@ namespace cheat{
// method to return a subset of allhits that are matched to a list of TrackIDs
const std::vector<std::vector<art::Ptr<recob::OpHit>>> TrackIDsToOpHits(std::vector<art::Ptr<recob::OpHit>> const& allhits,
- std::vector<int> const& tkIDs);
+ std::vector<int> const& tkIDs);
// method to return the EveIDs of particles contributing ionization
// electrons to the identified hit
@@ -99,43 +99,43 @@ namespace cheat{
std::vector<double> SimSDPsToXYZ(std::vector<sim::SDP> const& ides);
// method to return the XYZ position of the weighted average energy deposition for a given hit
- std::vector<double> OpHitToXYZ(art::Ptr<recob::OpHit> const& hit);
+ std::vector<double> OpHitToXYZ(art::Ptr<recob::OpHit> const& hit);
// method to return the XYZ position of a space point (unweighted average XYZ of component hits).
/* std::vector<double> SpacePointToXYZ(art::Ptr<recob::SpacePoint> const& spt,
- art::Event const& evt,
- std::string const& label);*/
+ art::Event const& evt,
+ std::string const& label);*/
// method to return the XYZ position of a space point (unweighted average XYZ of component hits).
/* std::vector<double> SpacePointHitsToXYZ(art::PtrVector<recob::Hit> const& hits);*/
// method to return the fraction of hits in a collection that come from the specified Geant4 track ids
double OpHitCollectionPurity(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& hits);
+ std::vector< art::Ptr<recob::OpHit> > const& hits);
// method to return the fraction of all hits in an event from a specific set of Geant4 track IDs that are
// represented in a collection of hits
double OpHitCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& hits,
- std::vector< art::Ptr<recob::OpHit> > const& allhits);
+ std::vector< art::Ptr<recob::OpHit> > const& hits,
+ std::vector< art::Ptr<recob::OpHit> > const& allhits);
double OpHitCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& hits,
- std::vector< art::Ptr<recob::OpHit> > const& allhits,
- geo::View_t const& view);
+ std::vector< art::Ptr<recob::OpHit> > const& hits,
+ std::vector< art::Ptr<recob::OpHit> > const& allhits,
+ geo::View_t const& view);
// method to return the fraction of charge in a collection that come from the specified Geant4 track ids
double OpHitChargeCollectionPurity(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& hits);
+ std::vector< art::Ptr<recob::OpHit> > const& hits);
// method to return the fraction of all charge in an event from a specific set of Geant4 track IDs that are
// represented in a collection of hits
double OpHitChargeCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& hits,
- std::vector< art::Ptr<recob::OpHit> > const& allhits);
+ std::vector< art::Ptr<recob::OpHit> > const& hits,
+ std::vector< art::Ptr<recob::OpHit> > const& allhits);
double OpHitChargeCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& hits,
- std::vector< art::Ptr<recob::OpHit> > const& allhits,
- geo::View_t const& view);
+ std::vector< art::Ptr<recob::OpHit> > const& hits,
+ std::vector< art::Ptr<recob::OpHit> > const& allhits,
+ geo::View_t const& view);
// method to return all EveIDs corresponding to the current sim::ParticleList
std::set<int> GetSetOfEveIDs();
@@ -152,9 +152,9 @@ namespace cheat{
const std::vector<const sim::OpDetBacktrackerRecord*>& OpDetBacktrackerRecords() const { return cOpDetBacktrackerRecords; }
void ChannelToTrackSDPs(std::vector<sim::TrackSDP>& trackSDPs,
- int channel,
- const double hit_start_time,
- const double hit_end_time);
+ int channel,
+ const double hit_start_time,
+ const double hit_end_time);
private:
larsim/MCCheater/PhotonBackTracker_service.cc
diff --git a/larsim/MCCheater/PhotonBackTracker_service.cc b/larsim/MCCheater/PhotonBackTracker_service.cc
index 6e073c8..b6b5b6d 100644
--- a/larsim/MCCheater/PhotonBackTracker_service.cc
+++ b/larsim/MCCheater/PhotonBackTracker_service.cc
@@ -31,7 +31,7 @@ namespace cheat{
//----------------------------------------------------------------------
PhotonBackTracker::PhotonBackTracker(const fhicl::ParameterSet& pset,
- art::ActivityRegistry& reg)
+ art::ActivityRegistry& reg)
{
reconfigure(pset);
@@ -67,8 +67,8 @@ namespace cheat{
// if that is the case, we'll take care of it later
if(pHandle.failedToGet()){
mf::LogWarning("PhotonBackTracker") << "failed to get handle to simb::MCParticle from "
- << fG4ModuleLabel
- << ", return";
+ << fG4ModuleLabel
+ << ", return";
return;
}
@@ -82,31 +82,31 @@ namespace cheat{
if( fo.isValid() ){
for(size_t p = 0; p < pHandle->size(); ++p){
-
- simb::MCParticle *part = new simb::MCParticle(pHandle->at(p));
- fParticleList.Add(part);
-
- // get the simb::MCTruth associated to this sim::ParticleList
- try{
- art::Ptr<simb::MCTruth> mct = fo.at(p);
- if(fMCTruthList.size() < 1) fMCTruthList.push_back(mct);
- else{
- // check that we are not adding a simb::MCTruth twice to the collection
- // we know that all the particles for a given simb::MCTruth are put into the
- // collection of particles at the same time, so we can just check that the
- // current art::Ptr has a different id than the last one put
- if(!(mct == fMCTruthList.back())) fMCTruthList.push_back(mct);
- }
- // fill the track id to mctruth index map
- fTrackIDToMCTruthIndex[pHandle->at(p).TrackId()] = fMCTruthList.size() - 1;
- }
- catch(cet::exception &ex){
- mf::LogWarning("PhotonBackTracker") << "unable to find MCTruth from ParticleList "
- << "created in " << fG4ModuleLabel << " "
- << "any attempt to get the MCTruth objects from "
- << "the backtracker will fail\n"
- << "message from caught exception:\n" << ex;
- }
+
+ simb::MCParticle *part = new simb::MCParticle(pHandle->at(p));
+ fParticleList.Add(part);
+
+ // get the simb::MCTruth associated to this sim::ParticleList
+ try{
+ art::Ptr<simb::MCTruth> mct = fo.at(p);
+ if(fMCTruthList.size() < 1) fMCTruthList.push_back(mct);
+ else{
+ // check that we are not adding a simb::MCTruth twice to the collection
+ // we know that all the particles for a given simb::MCTruth are put into the
+ // collection of particles at the same time, so we can just check that the
+ // current art::Ptr has a different id than the last one put
+ if(!(mct == fMCTruthList.back())) fMCTruthList.push_back(mct);
+ }
+ // fill the track id to mctruth index map
+ fTrackIDToMCTruthIndex[pHandle->at(p).TrackId()] = fMCTruthList.size() - 1;
+ }
+ catch(cet::exception &ex){
+ mf::LogWarning("PhotonBackTracker") << "unable to find MCTruth from ParticleList "
+ << "created in " << fG4ModuleLabel << " "
+ << "any attempt to get the MCTruth objects from "
+ << "the backtracker will fail\n"
+ << "message from caught exception:\n" << ex;
+ }
}// end loop over particles to get MCTruthList
}// end if fo.isValid()
@@ -119,10 +119,10 @@ namespace cheat{
fParticleList.AdoptEveIdCalculator(new sim::EmEveIdCalculator);
LOG_DEBUG("PhotonBackTracker") << "PhotonBackTracker has " << cOpDetBacktrackerRecords.size()
- << " sim::OpDetBacktrackerRecords and " << GetSetOfTrackIDs().size()
- << " tracks. The particles are:\n"
- << fParticleList
- << "\n the MCTruth information is\n";
+ << " sim::OpDetBacktrackerRecords and " << GetSetOfTrackIDs().size()
+ << " tracks. The particles are:\n"
+ << fParticleList
+ << "\n the MCTruth information is\n";
for(size_t mc = 0; mc < fMCTruthList.size(); ++mc)
LOG_DEBUG("PhotonBackTracker") << *(fMCTruthList.at(mc).get());
@@ -136,8 +136,8 @@ namespace cheat{
if(part_it == fParticleList.end()){
mf::LogWarning("PhotonBackTracker") << "can't find particle with track id "
- << id << " in sim::ParticleList"
- << " returning null pointer";
+ << id << " in sim::ParticleList"
+ << " returning null pointer";
return 0;
}
@@ -161,8 +161,8 @@ namespace cheat{
if(/* mct < 0 || */ mct > fMCTruthList.size() )
throw cet::exception("PhotonBackTracker") << "attempting to find MCTruth index for "
- << "out of range value: " << mct
- << "/" << fMCTruthList.size() << "\n";
+ << "out of range value: " << mct
+ << "/" << fMCTruthList.size() << "\n";
return fMCTruthList[mct];
}
@@ -179,12 +179,12 @@ namespace cheat{
// loop over the SDPMAP
for(auto mapitr = pdTimeSDPmap.begin(); mapitr != pdTimeSDPmap.end(); mapitr++){
-
- // loop over the vector of SDP objects.
- const std::vector<sim::SDP>& sdpvec = (*mapitr).second;
- for(size_t iv = 0; iv < sdpvec.size(); ++iv){
- if( abs(sdpvec[iv].trackID) == id) sdps.push_back(sdpvec[iv]);
- }
+
+ // loop over the vector of SDP objects.
+ const std::vector<sim::SDP>& sdpvec = (*mapitr).second;
+ for(size_t iv = 0; iv < sdpvec.size(); ++iv){
+ if( abs(sdpvec[iv].trackID) == id) sdps.push_back(sdpvec[iv]);
+ }
} // end loop over map from sim::OpDetBacktrackerRecord
} // end loop over sim::OpDetBacktrackerRecords
@@ -221,14 +221,14 @@ namespace cheat{
const double start = (pTime-pWidth)*1000-fDelay;
const double end = (pTime+pWidth)*1000-fDelay;
- this->ChannelToTrackSDPs(trackSDPs, opHit->OpChannel(), start, end);
+ this->ChannelToTrackSDPs(trackSDPs, opHit->OpChannel(), start, end);
return trackSDPs;
}
//----------------------------------------------------------------------
const std::vector<std::vector<art::Ptr<recob::OpHit>>> PhotonBackTracker::TrackIDsToOpHits(std::vector<art::Ptr<recob::OpHit>> const& allOpHits,
- std::vector<int> const& tkIDs)
+ std::vector<int> const& tkIDs)
{
// returns a subset of the opHits in the allOpHits collection that are matched
// to MC particles listed in tkIDs
@@ -369,7 +369,7 @@ namespace cheat{
// loop over the ides and extract the track ids
for(size_t i = 0; i < trackSDPs.size(); ++i) {
- trackIDs.insert(trackSDPs[i].trackID);
+ trackIDs.insert(trackSDPs[i].trackID);
}
itr++;
@@ -380,7 +380,7 @@ namespace cheat{
//----------------------------------------------------------------------
double PhotonBackTracker::OpHitCollectionPurity(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& opHits)
+ std::vector< art::Ptr<recob::OpHit> > const& opHits)
{
// get the list of EveIDs that correspond to the opHits in this collection
// if the EveID shows up in the input list of trackIDs, then it counts
@@ -397,10 +397,10 @@ namespace cheat{
// don't double count if this opHit has more than one of the
// desired track IDs associated with it
for(size_t e = 0; e < opHitTrackIDs.size(); ++e){
- if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end()){
- desired += 1.;
- break;
- }
+ if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end()){
+ desired += 1.;
+ break;
+ }
}
}// end loop over opHits
@@ -413,7 +413,7 @@ namespace cheat{
//----------------------------------------------------------------------
double PhotonBackTracker::OpHitChargeCollectionPurity(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& opHits)
+ std::vector< art::Ptr<recob::OpHit> > const& opHits)
{
// get the list of EveIDs that correspond to the opHits in this collection
// if the EveID shows up in the input list of trackIDs, then it counts
@@ -432,10 +432,10 @@ namespace cheat{
// don't double count if this opHit has more than one of the
// desired track IDs associated with it
for(size_t e = 0; e < opHitTrackIDs.size(); ++e){
- if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end()){
- desired += opHit->Area();
- break;
- }
+ if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end()){
+ desired += opHit->Area();
+ break;
+ }
}
}// end loop over opHits
@@ -449,16 +449,16 @@ namespace cheat{
//----------------------------------------------------------------------
double PhotonBackTracker::OpHitCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& opHits,
- std::vector< art::Ptr<recob::OpHit> > const& allOpHits,
- geo::View_t const& view)
+ std::vector< art::Ptr<recob::OpHit> > const& opHits,
+ std::vector< art::Ptr<recob::OpHit> > const& allOpHits,
+ geo::View_t const& view)
{
throw cet::exception("PhotonBackTracker")<<"This function is not supported. OpHits do not have type View.\n";
}
double PhotonBackTracker::OpHitCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& opHits,
- std::vector< art::Ptr<recob::OpHit> > const& allOpHits)
+ std::vector< art::Ptr<recob::OpHit> > const& opHits,
+ std::vector< art::Ptr<recob::OpHit> > const& allOpHits)
{
// get the list of EveIDs that correspond to the opHits in this collection
// and the energy associated with the desired trackID
@@ -478,11 +478,11 @@ namespace cheat{
// also don't double count if this opHit has more than one of the
// desired track IDs associated with it
for(size_t e = 0; e < opHitTrackIDs.size(); ++e){
- if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
- opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
- desired += 1.;
- break;
- }
+ if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
+ opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
+ desired += 1.;
+ break;
+ }
}
}// end loop over opHits
@@ -498,15 +498,15 @@ namespace cheat{
std::vector<sim::TrackSDP> opHitTrackIDs = this->OpHitToTrackID(opHit);
for(size_t e = 0; e < opHitTrackIDs.size(); ++e){
- // don't worry about opHits where the energy fraction for the chosen
- // trackID is < 0.1
- // also don't double count if this opHit has more than one of the
- // desired track IDs associated with it
- if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
- opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
- total += 1.;
- break;
- }
+ // don't worry about opHits where the energy fraction for the chosen
+ // trackID is < 0.1
+ // also don't double count if this opHit has more than one of the
+ // desired track IDs associated with it
+ if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
+ opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
+ total += 1.;
+ break;
+ }
}
}// end loop over all opHits
@@ -519,15 +519,15 @@ namespace cheat{
//----------------------------------------------------------------------
double PhotonBackTracker::OpHitChargeCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& opHits,
- std::vector< art::Ptr<recob::OpHit> > const& allOpHits,
- geo::View_t const& view)
+ std::vector< art::Ptr<recob::OpHit> > const& opHits,
+ std::vector< art::Ptr<recob::OpHit> > const& allOpHits,
+ geo::View_t const& view)
{
throw cet::exception("PhotonBackTracker")<<"This function is not supported. OpHits do not have type View.\n";
}
double PhotonBackTracker::OpHitChargeCollectionEfficiency(std::set<int> trackIDs,
- std::vector< art::Ptr<recob::OpHit> > const& opHits,
- std::vector< art::Ptr<recob::OpHit> > const& allOpHits)
+ std::vector< art::Ptr<recob::OpHit> > const& opHits,
+ std::vector< art::Ptr<recob::OpHit> > const& allOpHits)
{
// get the list of EveIDs that correspond to the opHits in this collection
// and the energy associated with the desired trackID
@@ -547,11 +547,11 @@ namespace cheat{
// also don't double count if this opHit has more than one of the
// desired track IDs associated with it
for(size_t e = 0; e < opHitTrackIDs.size(); ++e){
- if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
- opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
- desired += opHit->Area();
- break;
- }
+ if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
+ opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
+ desired += opHit->Area();
+ break;
+ }
}
}// end loop over opHits
@@ -567,15 +567,15 @@ namespace cheat{
std::vector<sim::TrackSDP> opHitTrackIDs = this->OpHitToTrackID(opHit);
for(size_t e = 0; e < opHitTrackIDs.size(); ++e){
- // don't worry about opHits where the energy fraction for the chosen
- // trackID is < 0.1
- // also don't double count if this opHit has more than one of the
- // desired track IDs associated with it
- if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
- opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
- total += opHit->Area();
- break;
- }
+ // don't worry about opHits where the energy fraction for the chosen
+ // trackID is < 0.1
+ // also don't double count if this opHit has more than one of the
+ // desired track IDs associated with it
+ if(trackIDs.find(opHitTrackIDs[e].trackID) != trackIDs.end() &&
+ opHitTrackIDs[e].energyFrac >= fMinOpHitEnergyFraction){
+ total += opHit->Area();
+ break;
+ }
}
}// end loop over all opHits
@@ -599,16 +599,16 @@ namespace cheat{
if(!opDet)
throw cet::exception("PhotonBackTracker") << "No sim::OpDetBacktrackerRecord corresponding "
- << "to opDetNum: " << opDetNum << "\n";
+ << "to opDetNum: " << opDetNum << "\n";
return opDet;
}
//----------------------------------------------------------------------
void PhotonBackTracker::ChannelToTrackSDPs(std::vector<sim::TrackSDP>& trackSDPs,
- int channel,
- const double opHit_start_time,
- const double opHit_end_time)
+ int channel,
+ const double opHit_start_time,
+ const double opHit_end_time)
{
trackSDPs.clear();
@@ -629,7 +629,7 @@ namespace cheat{
// first get the total energy represented by all track ids for
// this channel and range of tdc values
for(size_t e = 0; e < simSDPs.size(); ++e)
- totalE += simSDPs[e].energy;
+ totalE += simSDPs[e].energy;
// protect against a divide by zero below
@@ -638,21 +638,21 @@ namespace cheat{
// loop over the entries in the map and fill the input vectors
for(size_t e = 0; e < simSDPs.size(); ++e){
-
- if(simSDPs[e].trackID == sim::NoParticleId) continue;
-
- sim::TrackSDP info;
- info.trackID = simSDPs[e].trackID;
- info.energyFrac = simSDPs[e].energy/totalE;
- info.energy = simSDPs[e].energy;
-
- trackSDPs.push_back(info);
-
+
+ if(simSDPs[e].trackID == sim::NoParticleId) continue;
+
+ sim::TrackSDP info;
+ info.trackID = simSDPs[e].trackID;
+ info.energyFrac = simSDPs[e].energy/totalE;
+ info.energy = simSDPs[e].energy;
+
+ trackSDPs.push_back(info);
+
}
}// end try
catch(cet::exception e){
mf::LogWarning("PhotonBackTracker") << "caught exception \n"
- << e;
+ << e;
}
return;
@@ -696,12 +696,12 @@ namespace cheat{
z += weight * sdp.z;
}// end loop over sim::SDPs
-
+
// if the sum of the weights is still 0, then return
// the obviously stupid default values
if(w < 1.e-5)
throw cet::exception("PhotonBackTracker") << "No sim::SDPs providing non-zero number of photons"
- << " can't determine originating location from truth\n";
+ << " can't determine originating location from truth\n";
xyz[0] = x/w;
xyz[1] = y/w;
larsim/MCCheater/photonbacktracker.fcl
diff --git a/larsim/MCCheater/photonbacktracker.fcl b/larsim/MCCheater/photonbacktracker.fcl
index b02c251..3934c9a 100644
--- a/larsim/MCCheater/photonbacktracker.fcl
+++ b/larsim/MCCheater/photonbacktracker.fcl
@@ -5,13 +5,7 @@ standard_photonbacktracker:
G4ModuleLabel: "largeant" # module that produced the sim::Particle and sim::SimChannel objects
MinimumHitEnergyFraction: 0.1 # minimum fraction of energy a G4 trackID contributes to a hit to be
# counted in hit based efficiency and purity calculations
+ Delay: 0 #This number is the difference between when light arrives at the detector, and the time recorded in OpHits. This number is experiment specific and should be set by each experiment.
}
-jp250L_photonbacktracker: @local::standard_photonbacktracker
-bo_photonbacktracker: @local::standard_photonbacktracker
-argoneut_photonbacktracker: @local::standard_photonbacktracker
-microboone_photonbacktracker: @local::standard_photonbacktracker
-dune35t_photonbacktracker: @local::standard_photonbacktracker
-dunefd_photonbacktracker: @local::standard_photonbacktracker
-
END_PROLOG
larreco
larreco/HitFinder/GausHitFinder_module.cc
diff --git a/larreco/HitFinder/GausHitFinder_module.cc b/larreco/HitFinder/GausHitFinder_module.cc
index b35c72a..3273bd0 100644
--- a/larreco/HitFinder/GausHitFinder_module.cc
+++ b/larreco/HitFinder/GausHitFinder_module.cc
@@ -125,7 +125,7 @@ namespace hit{
size_t fMaxMultiHit; ///<maximum hits for multi fit
int fAreaMethod; ///<Type of area calculation
std::vector<double> fAreaNorms; ///<factors for converting area to same units as peak height
- bool fTryNplus1Fits; ///<whether we will (0) or won't (1) try n+1 fits
+ bool fTryNplus1Fits; ///<whether we will (true) or won't (false) try n+1 fits
double fChi2NDFRetry; ///<Value at which a second n+1 Fit will be tried
double fChi2NDF; ///maximum Chisquared / NDF allowed for a hit to be saved
size_t fNumBinsToAverage; ///< If bin averaging for peak finding, number bins to average
larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.cxx
diff --git a/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.cxx b/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.cxx
index ed5a357..15c78f8 100644
--- a/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.cxx
+++ b/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.cxx
@@ -235,15 +235,17 @@ bool nnet::DataProviderAlg::bufferPatch(size_t wire, float drift) const
int d0 = fCurrentScaledDrift - halfSizeD;
int d1 = fCurrentScaledDrift + halfSizeD;
+ int wsize = fWireDriftData.size();
for (int w = w0, wpatch = 0; w < w1; ++w, ++wpatch)
{
auto & dst = fWireDriftPatch[wpatch];
- if ((w >= 0) && (w < (int)fWireDriftData.size()))
+ if ((w >= 0) && (w < wsize))
{
auto & src = fWireDriftData[w];
+ int dsize = src.size();
for (int d = d0, dpatch = 0; d < d1; ++d, ++dpatch)
{
- if ((d >= 0) && (d < (int)src.size()))
+ if ((d >= 0) && (d < dsize))
{
dst[dpatch] = src[d];
}
@@ -299,12 +301,30 @@ bool nnet::DataProviderAlg::isInsideFiducialRegion(unsigned int wire, float drif
}
// ------------------------------------------------------
+
+
+// ------------------------------------------------------
+// -------------------ModelInterface---------------------
+// ------------------------------------------------------
+
+std::string nnet::ModelInterface::findFile(const char* fileName) const
+{
+ std::string fname_out;
+ cet::search_path sp("FW_SEARCH_PATH");
+ if (!sp.find_file(fileName, fname_out))
+ {
+ throw art::Exception(art::errors::NotFound)
+ << "Could not find the model file " << fileName;
+ }
+ return fname_out;
+}
+
// ------------------------------------------------------
// -----------------MlpModelInterface--------------------
// ------------------------------------------------------
nnet::MlpModelInterface::MlpModelInterface(const char* xmlFileName) :
- m(xmlFileName)
+ m(nnet::ModelInterface::findFile(xmlFileName).c_str())
{
mf::LogInfo("MlpModelInterface") << "MLP model loaded.";
}
@@ -357,7 +377,7 @@ float nnet::MlpModelInterface::GetOneOutput(int neuronIndex) const
// ------------------------------------------------------
nnet::KerasModelInterface::KerasModelInterface(const char* modelFileName) :
- m(modelFileName)
+ m(nnet::ModelInterface::findFile(modelFileName).c_str())
{
mf::LogInfo("KerasModelInterface") << "Keras model loaded.";
}
@@ -649,6 +669,7 @@ bool nnet::TrainingDataAlg::setWireEdepsAndLabels(
int best_pdg = pdgs[i0] & nnet::TrainingDataAlg::kPdgMask;
int vtx_flags = pdgs[i0] & nnet::TrainingDataAlg::kVtxMask;
+ int type_flags = pdgs[i0] & nnet::TrainingDataAlg::kTypeMask;
float max_edep = edeps[i0];
for (size_t k = i0 + 1; k < i1; ++k)
{
@@ -658,11 +679,13 @@ bool nnet::TrainingDataAlg::setWireEdepsAndLabels(
max_edep = ek;
best_pdg = pdgs[k] & nnet::TrainingDataAlg::kPdgMask; // remember best matching pdg
}
+ type_flags |= pdgs[k] & nnet::TrainingDataAlg::kTypeMask; // accumulate track type flags
vtx_flags |= pdgs[k] & nnet::TrainingDataAlg::kVtxMask; // accumulate all vtx flags
}
wEdep[i] = max_edep;
+ best_pdg |= type_flags;
if (fSaveVtxFlags) best_pdg |= vtx_flags;
wPdg[i] = best_pdg;
}
@@ -701,6 +724,32 @@ nnet::TrainingDataAlg::WireDrift nnet::TrainingDataAlg::getProjection(double x,
}
// ------------------------------------------------------
+bool nnet::TrainingDataAlg::isMuonDecaying(const simb::MCParticle & particle,
+ const std::unordered_map< int, const simb::MCParticle* > & particleMap) const
+{
+ bool hasElectron = false, hasNuMu = false, hasNuE = false;
+
+ int pdg = abs(particle.PdgCode());
+ if ((pdg == 13) && (particle.EndProcess() == "FastScintillation")) // potential muon decay at rest
+ {
+ unsigned int nSec = particle.NumberDaughters();
+ for (size_t d = 0; d < nSec; ++d)
+ {
+ auto d_search = particleMap.find(particle.Daughter(d));
+ if (d_search != particleMap.end())
+ {
+ auto const & daughter = *((*d_search).second);
+ int d_pdg = abs(daughter.PdgCode());
+ if (d_pdg == 11) hasElectron = true;
+ else if (d_pdg == 14) hasNuMu = true;
+ else if (d_pdg == 12) hasNuE = true;
+ }
+ }
+ }
+
+ return (hasElectron && hasNuMu && hasNuE);
+}
+
void nnet::TrainingDataAlg::collectVtxFlags(
std::unordered_map< size_t, std::unordered_map< int, int > > & wireToDriftToVtxFlags,
const std::unordered_map< int, const simb::MCParticle* > & particleMap,
@@ -729,26 +778,11 @@ void nnet::TrainingDataAlg::collectVtxFlags(
break;
case 13: // mu+/-
- if ((particle.EndProcess() == "FastScintillation")) // potential decay at rest
- {
- unsigned int nSec = particle.NumberDaughters();
-
- for (size_t d = 0; d < nSec; ++d)
- {
- auto d_search = particleMap.find(particle.Daughter(d));
- if (d_search != particleMap.end())
- {
- auto const & daughter = *((*d_search).second);
- int d_pdg = abs(daughter.PdgCode());
- if (d_pdg == 11)
- {
- //std::cout << "---> mu decay to electron" << std::endl;
- flagsEnd = nnet::TrainingDataAlg::kDecay;
- break;
- }
- }
- }
- }
+ if (nnet::TrainingDataAlg::isMuonDecaying(particle, particleMap))
+ {
+ //std::cout << "---> mu decay to electron" << std::endl;
+ flagsEnd = nnet::TrainingDataAlg::kDecay;
+ }
break;
case 111: // pi0
@@ -941,19 +975,49 @@ bool nnet::TrainingDataAlg::setEventData(const art::Event& event,
{
int pdg = 0;
int tid = energyDeposit.trackID;
- if (tid < 0) { pdg = 11; tid = -tid; } // negative tid means it is EM activity
-
- auto search = particleMap.find(tid);
- if (search == particleMap.end())
+ if (tid < 0) // negative tid means it is EM activity, and -tid is the mother
{
- mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
- continue;
+ pdg = 11; tid = -tid;
+
+ auto search = particleMap.find(tid);
+ if (search == particleMap.end())
+ {
+ mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
+ continue;
+ }
+ auto const & mother = *((*search).second); // mother particle of this EM
+ int mPdg = abs(mother.PdgCode());
+ if ((mPdg == 13) || (mPdg == 211) || (mPdg == 2212))
+ {
+ if (energyDeposit.numElectrons > 10) pdg |= nnet::TrainingDataAlg::kDelta; // tag delta ray
+ }
+ }
+ else
+ {
+ auto search = particleMap.find(tid);
+ if (search == particleMap.end())
+ {
+ mf::LogWarning("TrainingDataAlg") << "PARTICLE NOT FOUND";
+ continue;
+ }
+ auto const & particle = *((*search).second);
+ pdg = abs(particle.PdgCode());
+
+ if (pdg == 11) // electron, check if it is Michel
+ {
+ auto msearch = particleMap.find(particle.Mother());
+ if (msearch != particleMap.end())
+ {
+ auto const & mother = *((*msearch).second);
+ if (nnet::TrainingDataAlg::isMuonDecaying(mother, particleMap))
+ {
+ pdg |= nnet::TrainingDataAlg::kMichel; // tag Michel
+ }
+ }
+ }
}
- auto const & particle = *((*search).second);
- if (!pdg) pdg = particle.PdgCode(); // not EM activity so read what PDG it is
-
- trackToPDG[energyDeposit.trackID] = abs(pdg);
+ trackToPDG[energyDeposit.trackID] = pdg;
double energy = energyDeposit.numElectrons * electronsToGeV;
timeToTrackToCharge[time][energyDeposit.trackID] += energy;
@@ -962,6 +1026,7 @@ bool nnet::TrainingDataAlg::setEventData(const art::Event& event,
} // loop over time slices
} // for each SimChannel
+ int type_pdg_mask = nnet::TrainingDataAlg::kTypeMask | nnet::TrainingDataAlg::kPdgMask;
for (auto const & ttc : timeToTrackToCharge)
{
float max_deposit = 0.0;
@@ -974,11 +1039,11 @@ bool nnet::TrainingDataAlg::setEventData(const art::Event& event,
max_pdg = trackToPDG[tc.first];
}
}
-
+
if (ttc.first < (int)labels_deposit.size())
{
labels_deposit[ttc.first] = max_deposit;
- labels_pdg[ttc.first] = max_pdg & 0xFFFF;
+ labels_pdg[ttc.first] = max_pdg & type_pdg_mask;
}
}
larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.h
diff --git a/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.h b/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.h
index 109af89..720fe78 100644
--- a/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.h
+++ b/larreco/RecoAlg/ImagePatternAlgs/PointIdAlg/PointIdAlg.h
@@ -168,7 +168,7 @@ public:
protected:
ModelInterface(void) { }
-private:
+ std::string findFile(const char* fileName) const;
};
// ------------------------------------------------------
@@ -258,18 +258,28 @@ class nnet::TrainingDataAlg : public nnet::DataProviderAlg
{
public:
+ enum EMask
+ {
+ kNone = 0,
+ kPdgMask = 0x00000FFF, // pdg code mask
+ kTypeMask = 0x0000F000, // track type mask
+ kVtxMask = 0x0FFF0000 // vertex flags, still can use 0xFFFF0000 if more vtx types needed
+ };
+
+ enum ETrkType
+ {
+ kDelta = 0x1000, // delta electron
+ kMichel = 0x2000 // Michel electron
+ };
+
enum EVtxId
{
- kNone = 0,
kNuNC = 0x0010000, kNuCC = 0x0020000, // nu interaction type
kNuE = 0x0100000, kNuMu = 0x0200000, kNuTau = 0x0400000, // nu flavor
kHadr = 0x1000000, // hadronic inelastic scattering
kPi0 = 0x2000000, // pi0 produced in this vertex
kDecay = 0x4000000, // point of particle decay
kConv = 0x8000000, // gamma conversion
-
- kVtxMask = 0x0FFF0000, // vertex flags, still can use 0xFFFF0000 if more vtx types needed
- kPdgMask = 0x0000FFFF // pdg code mask (pdg is saved in the same int as vtx flags)
};
struct Config : public nnet::DataProviderAlg::Config
@@ -336,6 +346,10 @@ private:
const std::unordered_map< int, const simb::MCParticle* > & particleMap,
unsigned int view) const;
+ bool isMuonDecaying(
+ const simb::MCParticle & particle,
+ const std::unordered_map< int, const simb::MCParticle* > & particleMap) const;
+
std::vector< std::vector<float> > fWireDriftEdep;
std::vector< std::vector<int> > fWireDriftPdg;
larreco/RecoAlg/ImagePatternAlgs/scripts/utils.py
diff --git a/larreco/RecoAlg/ImagePatternAlgs/scripts/utils.py b/larreco/RecoAlg/ImagePatternAlgs/scripts/utils.py
index e620236..c8229a1 100644
--- a/larreco/RecoAlg/ImagePatternAlgs/scripts/utils.py
+++ b/larreco/RecoAlg/ImagePatternAlgs/scripts/utils.py
@@ -40,9 +40,9 @@ def get_data(fname, drift_margin = 0, crop = True):
A_pdg[deposit_th_ind] = 0
tracks = A_pdg.copy()
showers = A_pdg.copy()
- tracks[(A_pdg & 0xFFFF) == 11] = 0
+ tracks[(A_pdg & 0x0FFF) == 11] = 0
tracks[tracks > 0] = 1
- showers[(A_pdg & 0xFFFF) != 11] = 0
+ showers[(A_pdg & 0x0FFF) != 11] = 0
showers[showers > 0] = 1
return A_raw, A_deposit, A_pdg, tracks, showers
larreco/RecoAlg/TrackKalmanFitter.cxx
diff --git a/larreco/RecoAlg/TrackKalmanFitter.cxx b/larreco/RecoAlg/TrackKalmanFitter.cxx
index 19fe252..5bb6758 100644
--- a/larreco/RecoAlg/TrackKalmanFitter.cxx
+++ b/larreco/RecoAlg/TrackKalmanFitter.cxx
@@ -201,7 +201,6 @@ bool trkf::TrackKalmanFitter::fitTrack(const recob::Track& track, const std::vec
for (auto itertrack = fittedTrack.getTrackMap().rbegin(); itertrack != fittedTrack.getTrackMap().rend(); ++itertrack) {
trkf::KHitTrack& fwdTrack = itertrack->second;
trkf::KHitWireX khit(dynamic_cast<const trkf::KHitWireX&>(*fwdTrack.getHit().get()));//need a non const copy in case we want to modify the error
- if (useRMS_) khit.setMeasError(khit.getMeasError()*khit.getHit()->RMS()*khit.getHit()->RMS()/(std::max(khit.getHit()->SigmaPeakTime()*khit.getHit()->SigmaPeakTime(),0.08333333333f)));
boost::optional<double> pdist = prop_->noise_prop(trf,khit.getMeasSurface(),trkf::Propagator::BACKWARD,true);
if (!pdist) {
larreco/TrackFinder/KalmanFilterFinalTrackFitter_module.cc
diff --git a/larreco/TrackFinder/KalmanFilterFinalTrackFitter_module.cc b/larreco/TrackFinder/KalmanFilterFinalTrackFitter_module.cc
index d6728d2..3bf23a1 100644
--- a/larreco/TrackFinder/KalmanFilterFinalTrackFitter_module.cc
+++ b/larreco/TrackFinder/KalmanFilterFinalTrackFitter_module.cc
@@ -197,6 +197,9 @@ trkf::KalmanFilterFinalTrackFitter::KalmanFilterFinalTrackFitter(trkf::KalmanFil
produces<std::vector<recob::Track> >();
produces<art::Assns<recob::Track, recob::Hit> >();
+ if (p_().options().trackFromPF()) {
+ produces<art::Assns<recob::PFParticle, recob::Track> >();
+ }
//throw expections to avoid possible silent failures due to incompatible configuration options
if (p_().options().trackFromPF()==0 && p_().options().idFromPF())
@@ -209,7 +212,7 @@ trkf::KalmanFilterFinalTrackFitter::KalmanFilterFinalTrackFitter(trkf::KalmanFil
throw cet::exception("KalmanFilterFinalTrackFitter") << "Incompatible configuration parameters: cannot use dirFromVec=true with dirFromVtxPF=true." << "\n";
unsigned int nPFroms = 0;
if (p_().options().pFromCalo()) nPFroms++;
- if (p_().options().pFromMSChi2()) nPFroms++;
+ if (p_().options().pFromMSChi2()) nPFroms++;
if (p_().options().pFromLength()) nPFroms++;
if (p_().options().pFromMC()) nPFroms++;
if (nPFroms>1) {
@@ -251,6 +254,8 @@ void trkf::KalmanFilterFinalTrackFitter::produce(art::Event & e)
}
if (p_().options().trackFromPF()) {
+
+ auto outputPFAssn = std::make_unique<art::Assns<recob::PFParticle, recob::Track> >();
art::ValidHandle<std::vector<recob::PFParticle> > inputPFParticle = e.getValidHandle<std::vector<recob::PFParticle> >(pfParticleInputTag);
assocTracks = std::unique_ptr<art::FindManyP<recob::Track> >(new art::FindManyP<recob::Track>(inputPFParticle, e, pfParticleInputTag));
@@ -297,8 +302,12 @@ void trkf::KalmanFilterFinalTrackFitter::produce(art::Event & e)
for (auto const& trhit: outHits) {
outputHits->addSingle(aptr, trhit);
}
+ outputPFAssn->addSingle(art::Ptr<recob::PFParticle>(inputPFParticle, iPF), aptr);
}
}
+ e.put(std::move(outputTracks));
+ e.put(std::move(outputHits));
+ e.put(std::move(outputPFAssn));
} else {
art::ValidHandle<std::vector<recob::Track> > inputTracks = e.getValidHandle<std::vector<recob::Track> >(trackInputTag);
@@ -342,11 +351,9 @@ void trkf::KalmanFilterFinalTrackFitter::produce(art::Event & e)
outputHits->addSingle(aptr, trhit);
}
}
+ e.put(std::move(outputTracks));
+ e.put(std::move(outputHits));
}
-
- e.put(std::move(outputTracks));
- e.put(std::move(outputHits));
-
}
double trkf::KalmanFilterFinalTrackFitter::setMomValue(art::Ptr<recob::Track> ptrack, const std::unique_ptr<art::FindManyP<anab::Calorimetry> >& trackCalo, const double pMC, const int pId) const {
larana
larana/OpticalDetector/OpHitFinder/UtilFunc.cxx
diff --git a/larana/OpticalDetector/OpHitFinder/UtilFunc.cxx b/larana/OpticalDetector/OpHitFinder/UtilFunc.cxx
index 8f6ec29..e09a8d8 100644
--- a/larana/OpticalDetector/OpHitFinder/UtilFunc.cxx
+++ b/larana/OpticalDetector/OpHitFinder/UtilFunc.cxx
@@ -66,7 +66,7 @@ namespace pmtana {
double bin_width = ((*res.second) - (*res.first)) / ((double)nbins);
- if(nbins==1) return ((*res.first) + bin_width /2.);
+ if(nbins==1 || bin_width == 0) return ((*res.first) + bin_width /2.);
//std::cout<<"Min: "<<(*res.first)<<" Max: "<<(*res.second)<<" Width: "<<bin_width<<std::endl;
larana/T0Finder/PhotonCounterT0Matching_module.cc
diff --git a/larana/T0Finder/PhotonCounterT0Matching_module.cc b/larana/T0Finder/PhotonCounterT0Matching_module.cc
index 9e1267a..4b4b4d8 100644
--- a/larana/T0Finder/PhotonCounterT0Matching_module.cc
+++ b/larana/T0Finder/PhotonCounterT0Matching_module.cc
@@ -125,6 +125,7 @@ private:
double fWeightOfDeltaYZ;
double fMatchCriteria;
double fPEThreshold;
+ bool fVerbosity;
// Variables used in module.......
std::vector<double> trackStart;
@@ -189,6 +190,8 @@ void lbne::PhotonCounterT0Matching::reconfigure(fhicl::ParameterSet const & p)
fWeightOfDeltaYZ = (p.get< double > ("WeightOfDeltaYZ" ) );
fMatchCriteria = (p.get< double > ("MatchCriteria" ) );
fPEThreshold = (p.get< double > ("PEThreshold" ) );
+
+ fVerbosity = (p.get< bool > ("Verbose",false) );
}
void lbne::PhotonCounterT0Matching::beginJob()
@@ -219,7 +222,7 @@ void lbne::PhotonCounterT0Matching::beginJob()
20, 0, 100, 60, 0, 300);
hFitParam_Length = tfs->make<TH2D>("hFitParam_Length", "How fit correlates with track length; Fit correlation; Track Length (cm)", 50, 0, 250, 30, 0, 300);
hPhotonT0_MCT0 = tfs->make<TH2D>("hPhotonT0_MCT0" , "Comparing Photon Counter reconstructed T0 against MCTruth T0; Photon Counter T0 (us); MCTruthT0 T0 (us)", 1760, -1600, 16000, 1760, -1600, 16000);
- hT0_diff_full = tfs->make<TH1D>("hT0_diff_full" , "Difference between MCTruth T0 and photon detector T0; Time difference (us); Number", 320, -1600, 1600);
+ hT0_diff_full = tfs->make<TH1D>("hT0_diff_full" , "Difference between MCTruth T0 and photon detector T0; Time difference (us); Number", 500, -2500, 2500);
hT0_diff_zoom = tfs->make<TH1D>("hT0_diff_zoom" , "Difference between MCTruth T0 and photon detector T0; Time difference (us); Number", 320, -1.6, 1.6);
}
@@ -269,11 +272,12 @@ void lbne::PhotonCounterT0Matching::produce(art::Event & evt)
size_t NTracks = tracklist.size();
size_t NFlashes = flashlist.size();
- std::cout << "There were " << NTracks << " tracks and " << NFlashes << " flashes in this event." << std::endl;
+ if (fVerbosity)
+ std::cout << "There were " << NTracks << " tracks and " << NFlashes << " flashes in this event." << std::endl;
// Now to access PhotonCounter for each track...
for(size_t iTrk=0; iTrk < NTracks; ++iTrk) {
- std::cout << "\n New Track " << (int)iTrk << std::endl;
+ if (fVerbosity) std::cout << "\n New Track " << (int)iTrk << std::endl;
// Reset Variables.
BestFlashTime = BestFitParam = BestTrackCentre_X = BestTrackLength = 9999;
BestTimeSepPredX = BestPredictedX = BestDeltaPredX = BestminYZSep = MCTruthT0 = 9999;
@@ -292,13 +296,13 @@ void lbne::PhotonCounterT0Matching::produce(art::Event & evt)
TrackLength);
// Some cout statement about track properties.
- ///*
- std::cout << trackStart[0] << " " << trackEnd[0] << " " << TrackLength_X << " " << TrackCentre_X
- << "\n" << trackStart[1] << " " << trackEnd[1] << " " << TrackLength_Y << " " << TrackCentre_Y
- << "\n" << trackStart[2] << " " << trackEnd[2] << " " << TrackLength_Z << " " << TrackCentre_Z
- << "\n" << trkTimeStart << " " << trkTimeEnd << " " << trkTimeLengh << " " << trkTimeCentre
- << std::endl;
- //*/
+ if (fVerbosity) {
+ std::cout << trackStart[0] << " " << trackEnd[0] << " " << TrackLength_X << " " << TrackCentre_X
+ << "\n" << trackStart[1] << " " << trackEnd[1] << " " << TrackLength_Y << " " << TrackCentre_Y
+ << "\n" << trackStart[2] << " " << trackEnd[2] << " " << TrackLength_Z << " " << TrackCentre_Z
+ << "\n" << trkTimeStart << " " << trkTimeEnd << " " << trkTimeLengh << " " << trkTimeCentre
+ << std::endl;
+ }
// ----- Loop over flashes ------
for ( size_t iFlash=0; iFlash < NFlashes; ++iFlash ) {
//Reset some flash specific quantities
@@ -334,8 +338,10 @@ void lbne::PhotonCounterT0Matching::produce(art::Event & evt)
else if (fMatchCriteria == 2) FitParam = DeltaPredX;
//----FLASH INFO-----
- std::cout << "\nFlash " << (int)iFlash << " " << TrackCentre_X << ", " << TimeSepPredX << " - " << PredictedX << " = " << DeltaPredX << ", " << minYZSep << " -> " << FitParam << std::endl;
-
+ if (fVerbosity) {
+ std::cout << "\nFlash " << (int)iFlash << " " << TrackCentre_X << ", " << TimeSepPredX << " - " << PredictedX << " = "
+ << DeltaPredX << ", " << minYZSep << " -> " << FitParam << std::endl;
+ }
//----Select best flash------
//double YFitRegion = (-1 * DeltaPredX ) + 80;
//if ( minYZSep > YFitRegion ) continue;