diff --git a/fcl/reco/Definitions/stage0_icarus_defs.fcl b/fcl/reco/Definitions/stage0_icarus_defs.fcl index 37841bca0..c22d2f9a9 100644 --- a/fcl/reco/Definitions/stage0_icarus_defs.fcl +++ b/fcl/reco/Definitions/stage0_icarus_defs.fcl @@ -9,6 +9,7 @@ #include "recowire_icarus.fcl" #include "hitfindermodules_icarus.fcl" #include "timing_icarus.fcl" +#include "timing_beam.fcl" #include "icarus_ophitfinder.fcl" #include "icarus_flashfinder.fcl" #include "trigger_emulation_icarus.fcl" @@ -115,6 +116,9 @@ icarus_stage0_producers: daqPMTonbeam: @local::copyPMTonBeam + ### Beam timing + beamTiming: @local::icarus_beam_signal_extractor + ### Purity monitoring purityana0: { module_type: "ICARUSPurityDQM" } purityana1: { module_type: "ICARUSPurityDQM" } @@ -243,6 +247,7 @@ icarus_stage0_PMT: [ triggerconfig, daqTrigger, pmtconfig, daqPMT, + beamTiming, pmtconfigbaselines, pmtthr, pmtbaselines, diff --git a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl index 41f8ea1b8..62ff01af7 100644 --- a/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl +++ b/fcl/reco/Stage1/Run2/stage1_run2_icarus_MC.fcl @@ -50,7 +50,7 @@ physics.producers: { physics.reco: [ @sequence::icarus_reco_Gauss_CryoE , @sequence::icarus_reco_Gauss_CryoW , - @sequence::icarus_reco_fm, + @sequence::icarus_reco_fm, @sequence::icarus_tpcpmtbarycentermatch, caloskimCalorimetryCryoE, caloskimCalorimetryCryoW, mcassociationsGausCryoE, mcassociationsGausCryoW @@ -82,6 +82,11 @@ physics.analyzers.caloskimW.RawDigitproducers: ["MCDecodeTPCROI:PHYSCRATEDATATPC physics.producers.mcassociationsGausCryoE.HitParticleAssociations.HitModuleLabelVec: ["cluster3DCryoE"] physics.producers.mcassociationsGausCryoW.HitParticleAssociations.HitModuleLabelVec: ["cluster3DCryoW"] +# Remove missing products in MC +physics.analyzers.simpleLightAna.TriggerLabel: "" +physics.analyzers.simpleLightAna.RWMLabel: "" +physics.analyzers.simpleLightAna.OpDetWaveformLabels: ["opdaq"] + services.message.destinations : { STDCOUT: diff --git a/icaruscode/IcarusObj/PMTBeamSignal.h b/icaruscode/IcarusObj/PMTBeamSignal.h new file mode 100644 index 000000000..842af0375 --- /dev/null +++ b/icaruscode/IcarusObj/PMTBeamSignal.h @@ -0,0 +1,69 @@ +/** + * @file icaruscode/IcarusObj/PMTBeamSignal.h + * @brief Holds the event-by-event RWM or EW times + * @author Matteo Vicenzi (mvicenzi@bnl.gov) + * @date March 14 2024 + */ + +#ifndef ICARUSCODE_ICARUSOBJ_PMTBEAMSIGNAL_H +#define ICARUSCODE_ICARUSOBJ_PMTBEAMSIGNAL_H + +// C/C++ standard libraries +#include +#include +#include + +namespace icarus::timing +{ + + /// Special value to denote no special channel information. + static constexpr auto NoChannel = std::numeric_limits::max(); + /// Special value to denote no time channel information. + static constexpr double NoTime = std::numeric_limits::max(); + // Special value to denote no sample information. + static constexpr std::size_t NoSample = 0; + + /** + * @brief Beam time as seen by a PMT readout board. + * + * This could either be an early warning (EW) or a resistive wall monitor (RWM) time. + * These signals are delivered via fibers and digitized in special PMT channels. + * + * Both the time in @ref DetectorClocksElectronicsTime "electronics time scale" + * and time time relative to the hardware trigger are included. + * + * The information in this object may be missing: its validity should + * always be checked in advance with `isValid()`. + */ + + struct PMTBeamSignal + { + + /// Special channel this time was extracted from. + /// These are defined in `CAEN_V1730_setup_icarus.fcl`. + unsigned int specialChannel = NoChannel; + + /// Board on which the special channel is on (e.g: WW-TOP-A). + /// Should match the same format as `icarusDB::PMTChannelInfo_t::digitizerLabel`. + std::string digitizerLabel = ""; + + /// Crate this time applies to (e.g.: WW-TOP). + /// Corresponds to the first part of `digitizerLabel`. + std::string crate = ""; + + /// Sample within the waveform where the reference signal is found. + std::size_t sample = NoSample; + + /// Start time in electronics time [us]. + double startTimeAbs = NoTime; + + /// Start time relative to trigger time [us]. + double startTime = NoTime; + + /// Returns whether the time is valid. + bool isValid() const { return (sample != NoSample); } + }; + +} // namespace icarus::timing + +#endif // ICARUSCODE_ICARUSOBJ_PMTBEAMSIGNAL_H diff --git a/icaruscode/PMT/OpReco/CMakeLists.txt b/icaruscode/PMT/OpReco/CMakeLists.txt index ac84f4f01..129cd515f 100644 --- a/icaruscode/PMT/OpReco/CMakeLists.txt +++ b/icaruscode/PMT/OpReco/CMakeLists.txt @@ -16,9 +16,10 @@ cet_build_plugin(ICARUSOpHitFinder art::module set( MODULE_LIBRARIES - icarusalg::Utilities + icarusalg::Utilities sbnobj::Common_Trigger - larcorealg::Geometry + icaruscode::Decode_DataProducts + larcorealg::Geometry larcore::Geometry_Geometry_service lardataobj::RecoBase lardataobj::Simulation @@ -48,7 +49,7 @@ set( MODULE_LIBRARIES ROOT::Gdml ROOT::FFTW ROOT::Core - ROOT::Tree + ROOT::Tree ) cet_build_plugin(FakeFlash art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(FakePhotoS art::module LIBRARIES ${MODULE_LIBRARIES}) @@ -61,6 +62,7 @@ cet_build_plugin(ICARUSOpHitAna art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(ICARUSOpHitTuple art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(ICARUSParticleAna art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(TPCPMTBarycenterMatchProducer art::module LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(ICARUSBeamStructureAna art::module LIBRARIES ${MODULE_LIBRARIES}) install_headers() diff --git a/icaruscode/PMT/OpReco/ICARUSBeamStructureAna_module.cc b/icaruscode/PMT/OpReco/ICARUSBeamStructureAna_module.cc new file mode 100644 index 000000000..6771e596c --- /dev/null +++ b/icaruscode/PMT/OpReco/ICARUSBeamStructureAna_module.cc @@ -0,0 +1,561 @@ +//////////////////////////////////////////////////////////////////////// +// Class: ICARUSBeamStructureAna +// Plugin Type: analyzer (Unknown Unknown) +// File: ICARUSBeamStructureAna_module.cc +// +// Generated at Mon Jul 22 13:50:03 2024 by Matteo Vicenzi using cetskelgen +// from version . +//////////////////////////////////////////////////////////////////////// + +// framework libraries +#include "canvas/Utilities/InputTag.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "canvas/Utilities/InputTag.h" +#include "canvas/Utilities/Exception.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "canvas/Persistency/Common/FindOne.h" +#include "canvas/Persistency/Common/Assns.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" +#include "art_root_io/TFileService.h" + +// LArSoft libraries +#include "icaruscode/Decode/DataProducts/TriggerConfiguration.h" +#include "lardataobj/RecoBase/OpFlash.h" +#include "lardataobj/RecoBase/OpHit.h" +#include "icaruscode/Decode/ChannelMapping/IICARUSChannelMap.h" +#include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom() +#include "lardataalg/DetectorInfo/DetectorTimingTypes.h" // electronics_time +#include "lardataobj/RawData/OpDetWaveform.h" +#include "lardataobj/Simulation/BeamGateInfo.h" +#include "lardataobj/RawData/TriggerData.h" +#include "sbnobj/Common/Trigger/ExtraTriggerInfo.h" +#include "icaruscode/CRT/CRTUtils/CRTPMTMatchingUtils.h" +#include "sbnobj/Common/CRT/CRTPMTMatching.hh" +#include "icaruscode/IcarusObj/PMTBeamSignal.h" + +// ROOT libraries +#include "TTree.h" +#include "TFile.h" + +// C/C++ standard libraries +#include +#include +#include +#include + +// ----------------------------------------------------------------------------- +namespace opana +{ + class ICARUSBeamStructureAna; +} + +class opana::ICARUSBeamStructureAna : public art::EDAnalyzer +{ +public: + struct Config + { + fhicl::Sequence FlashLabels{ + fhicl::Name("FlashLabels"), + fhicl::Comment("Tags for the recob::OpFlash data products")}; + fhicl::Atom TriggerLabel{ + fhicl::Name("TriggerLabel"), + fhicl::Comment("Tag for trigger info")}; + fhicl::Atom RWMLabel{ + fhicl::Name("RWMLabel"), + fhicl::Comment("Tag for RWM info")}; + fhicl::Atom TriggerConfigLabel{ + fhicl::Name("TriggerConfigLabel"), + fhicl::Comment("Trigger configuration label")}; + fhicl::Atom CRTPMTMatchingLabel{ + fhicl::Name("CRTPMTMatchingLabel"), + fhicl::Comment("CRTPMT matching label")}; + }; // struct Config + + using Parameters = art::EDAnalyzer::Table; + + /// constructor + explicit ICARUSBeamStructureAna(Parameters const &config); + + /// Return RWM-relative time from a trigger-relative time + double getRWMRelativeTime(int channel, double t); + + /// Return side/wall from channel id + int getSideByChannel(const int channel); + + /// Return the RWM-relative flash interaction time + double getFlashBunchTime(std::vector channels, std::vector hit_rise_time); + + /// Clear all data structures + void clear(); + + void analyze(art::Event const &e) override; + void beginJob(); + void beginRun(const art::Run &run) override; + +private: + art::ServiceHandle tfs; + icarus::TriggerConfiguration fTriggerConfiguration; + + std::vector fFlashLabels; + art::InputTag fTriggerLabel; + art::InputTag fRWMLabel; + art::InputTag fTriggerConfigurationLabel; + art::InputTag fCRTPMTMatchingLabel; + + /// data members + std::vector fOpFlashTrees; + + int m_run; + int m_event; + int m_timestamp; + + // trigger info + unsigned int m_gate_type; + int m_trigger_type = -1; + std::string m_gate_name; + uint64_t m_trigger_timestamp; + uint64_t m_beam_gate_timestamp; + double m_beam_us; + double m_trigger_us; + double m_beam_gate_width; + + // flash info + int m_cryo; + int m_flash_id; + double m_flash_time; + double m_flash_time_rwm; + double m_flash_z; + double m_flash_y; + double m_flash_pe; + int m_flash_nhits; + std::vector m_channel_id; + std::vector m_hit_start_time; + std::vector m_hit_peak_time; + std::vector m_hit_rise_time; + std::vector m_hit_start_time_rwm; + std::vector m_hit_peak_time_rwm; + std::vector m_hit_rise_time_rwm; + std::vector m_hit_pe; + + // crt-pmt match + int m_flash_classification; + int m_flash_ncrthits; + std::vector m_crthit_x; + std::vector m_crthit_y; + std::vector m_crthit_z; + std::vector m_crttime_us; + std::vector m_crtpmttimediff_ns; + std::vector m_crtsys; + std::vector m_crtregion; + + // RWM times + std::vector fRWMTimes; +}; + +// -------------------------------------------------------------------------- +opana::ICARUSBeamStructureAna::ICARUSBeamStructureAna(Parameters const &config) + : art::EDAnalyzer(config), + fFlashLabels(config().FlashLabels()), + fTriggerLabel(config().TriggerLabel()), + fRWMLabel(config().RWMLabel()), + fTriggerConfigurationLabel(config().TriggerConfigLabel()), + fCRTPMTMatchingLabel(config().CRTPMTMatchingLabel()) +{ +} + +// --------------------------------------------------------------------------- +void opana::ICARUSBeamStructureAna::beginJob() +{ + + if (!fFlashLabels.empty()) + { + + for (auto const &label : fFlashLabels) + { + + // TTree for the flash in a given cryostat + std::string name = "beamtiming_" + label.label(); + std::string info = "Beam timing from label " + label.label(); + + TTree *ttree = tfs->make(name.c_str(), info.c_str()); + ttree->Branch("run", &m_run); + ttree->Branch("event", &m_event); + ttree->Branch("timestamp", &m_timestamp); + + ttree->Branch("gate_type", &m_gate_type); + ttree->Branch("gate_name", &m_gate_name); + ttree->Branch("beam_gate_us", &m_beam_us); + ttree->Branch("trigger_us", &m_trigger_us); + ttree->Branch("beam_gate_width", &m_beam_gate_width); + ttree->Branch("trigger_type", &m_trigger_type, "trigger_type/I"); + ttree->Branch("trigger_timestamp", &m_trigger_timestamp, "trigger_timestamp/l"); + ttree->Branch("beam_gate_timestamp", &m_beam_gate_timestamp, "beam_gate_timestamp/l"); + + ttree->Branch("cryo", &m_cryo); + ttree->Branch("flash_id", &m_flash_id); + ttree->Branch("flash_time", &m_flash_time); + ttree->Branch("flash_time_rwm", &m_flash_time_rwm); + ttree->Branch("flash_pe", &m_flash_pe); + ttree->Branch("flash_z", &m_flash_z); + ttree->Branch("flash_y", &m_flash_y); + ttree->Branch("flash_nhits", &m_flash_nhits); + ttree->Branch("channels", &m_channel_id); + ttree->Branch("hit_start_time", &m_hit_start_time); + ttree->Branch("hit_peak_time", &m_hit_peak_time); + ttree->Branch("hit_rise_time", &m_hit_rise_time); + ttree->Branch("hit_start_time_rwm", &m_hit_start_time_rwm); + ttree->Branch("hit_peak_time_rwm", &m_hit_peak_time_rwm); + ttree->Branch("hit_rise_time_rwm", &m_hit_rise_time_rwm); + ttree->Branch("hit_pe", &m_hit_pe); + + ttree->Branch("flash_classification", &m_flash_classification); + ttree->Branch("flash_ncrthits", &m_flash_ncrthits); + ttree->Branch("crthit_x", &m_crthit_x); + ttree->Branch("crthit_y", &m_crthit_y); + ttree->Branch("crthit_z", &m_crthit_z); + ttree->Branch("crthit_sys", &m_crtsys); + ttree->Branch("crthit_region", &m_crtregion); + ttree->Branch("crttime_us", &m_crttime_us); + ttree->Branch("crtpmttimediff_ns", &m_crtpmttimediff_ns); + + fOpFlashTrees.push_back(ttree); + } + } + else + { + mf::LogError("ICARUSBeamStructureAna") + << "No flash labels selected!!"; + } +} + +// ------------------------------------------------------------------------------ + +void opana::ICARUSBeamStructureAna::beginRun(const art::Run &r) +{ + fTriggerConfiguration = r.getProduct(fTriggerConfigurationLabel); +} + +// ------------------------------------------------------------------------------- + +void opana::ICARUSBeamStructureAna::analyze(art::Event const &e) +{ + // ---- + // Event metadata information + m_run = e.id().run(); + m_event = e.id().event(); + m_timestamp = e.time().timeHigh(); // precision to the second + + // ---- + // Trigger metadata information + if (!fTriggerLabel.empty()) + { + + auto const &extraInfo = e.getProduct(fTriggerLabel); + + if (extraInfo.isValid()) + { + + sbn::triggerSource bit = extraInfo.sourceType; + m_gate_type = (unsigned int)bit; // 1 BNB 2 NumI 3 offbeamBNB 4 offbeamNuMi + m_gate_name = bitName(bit); + m_trigger_type = value(extraInfo.triggerType); // 1 majority, 2 minbias + + // absolute timestamp + m_trigger_timestamp = extraInfo.triggerTimestamp; + m_beam_gate_timestamp = extraInfo.beamGateTimestamp; + + // time in electronics time + m_trigger_us = e.getProduct>(fTriggerLabel).at(0).TriggerTime(); // us + m_beam_us = e.getProduct>(fTriggerLabel).at(0).BeamGateTime(); // us + m_beam_gate_width = fTriggerConfiguration.getGateWidth(m_gate_type); + } + else + { + mf::LogError("ICARUSBeamStructureAna") << "No raw::Trigger associated to label: " << fTriggerLabel.label() << "!"; + } + } + else + { + mf::LogError("ICARUSBeamStructureAna") << "No trigger labels selected!!"; + } + + // ---- + // RWM times + + fRWMTimes = e.getProduct>(fRWMLabel); + if (fRWMTimes.empty()) + mf::LogTrace("ICARUSBeamStructureAna") << "Data product std::vector for '" << fRWMLabel.label() + << "' is empty in " << m_gate_name << " event!"; + + // ---- + // FLASH/CRT timing information + + if (!fFlashLabels.empty()) + { + + for (std::size_t iFlashLabel = 0; iFlashLabel < fFlashLabels.size(); iFlashLabel++) + { + + auto const label = fFlashLabels[iFlashLabel]; + m_cryo = iFlashLabel; + + auto const &flash_handle = e.getValidHandle>(label); + auto const &flashes = *flash_handle; + + // we want our flashes to be valid and not empty + if (!flash_handle.isValid()) + { + mf::LogError("ICARUSBeamStructureAna") << "Not found a recob::OpFlash with label '" << label.encode() << "'"; + } + else if (flashes.empty()) + { + mf::LogWarning("ICARUSBeamStructureAna") << "No recob::OpFlash in collection with label '" << label.encode() << "'"; + } + else + { + + art::FindManyP ophitsPtr(flash_handle, e, label); + art::FindOneP matchPtr(flash_handle, e, fCRTPMTMatchingLabel); + + // loop all flashes + std::size_t idx = 0; + for (auto const &flash : flashes) + { + + // Filling flash info... + m_flash_id = idx; + m_flash_time = flash.Time(); + m_flash_pe = flash.TotalPE(); + m_flash_z = flash.ZCenter(); + m_flash_y = flash.YCenter(); + + // ---- + // CRT match info + + auto const &match = matchPtr.at(idx); + + // if there is no match, there is no product + // fill null parameters + if (!match) + { + m_flash_classification = 0; + m_flash_ncrthits = 0; + } + else + { + m_flash_classification = static_cast(match->flashClassification); + m_flash_ncrthits = match->matchedCRTHits.size(); + + for (auto const &crthit : match->matchedCRTHits) + { + m_crthit_x.push_back(crthit.position.X()); + m_crthit_y.push_back(crthit.position.Y()); + m_crthit_z.push_back(crthit.position.Z()); + m_crttime_us.push_back(crthit.time); + m_crtpmttimediff_ns.push_back(1e3 * crthit.PMTTimeDiff); + m_crtsys.push_back(crthit.sys); + m_crtregion.push_back(crthit.region); + } + } + + // ---- + // OPHITS info + + auto const &ophits = ophitsPtr.at(idx); + + std::map startmap; + std::map peakmap; + std::map risemap; + std::map pemap; + + // loop all hits in the flash: save only the first one + for (auto const hit : ophits) + { + + const int ch = hit->OpChannel(); + double ts = hit->StartTime(); + double tp = hit->PeakTime(); + double tr = hit->RiseTime(); + double pe = hit->PE(); + + // select the first ophit (by time) in each channel + if (startmap.find(ch) != startmap.end()) + { + if (ts < startmap[ch]) + { + startmap[ch] = ts; + peakmap[ch] = tp; + risemap[ch] = ts + tr; + pemap[ch] = pe; + } + } + else + { + startmap[ch] = ts; + peakmap[ch] = tp; + risemap[ch] = ts + tr; + pemap[ch] = pe; + } + } + + // get number of unique PMTs in flash + m_flash_nhits = startmap.size(); + + for (auto it = startmap.begin(); it != startmap.end(); it++) + { + int ch = it->first; + m_channel_id.push_back(ch); + m_hit_start_time.push_back(it->second); + m_hit_peak_time.push_back(peakmap[ch]); + m_hit_rise_time.push_back(risemap[ch]); + m_hit_start_time_rwm.push_back(getRWMRelativeTime(ch, it->second)); + m_hit_peak_time_rwm.push_back(getRWMRelativeTime(ch, peakmap[ch])); + m_hit_rise_time_rwm.push_back(getRWMRelativeTime(ch, risemap[ch])); + m_hit_pe.push_back(pemap[ch]); + } + + // get the flash interaction time w.r.t. RWM + // this is currently the mean between the first ophits on opposite walls + m_flash_time_rwm = getFlashBunchTime(m_channel_id, m_hit_rise_time_rwm); + + fOpFlashTrees[iFlashLabel]->Fill(); + + clear(); + idx++; + } + } + } + } +} + +// ----------------------------------------------------------------------------- + +int opana::ICARUSBeamStructureAna::getSideByChannel(const int channel) +{ + + /* + Channels are numbered from east to west, from North (cryo side) to South (beam side) + We look in the opposide direction wrt to the beam direction South->North: + - Left is the east wall of each cryostat; + - Right is the west side of each cryostat; + - [ 0:89 ] and [180:269] are on the left, + the return value of the function is 0; + - [ 90-179 ] and [ 270:359 ] are on the right, + the return value of the function is 1; + */ + + int side = channel / 90; // always round down + return side % 2; +} + +// ----------------------------------------------------------------------------- + +double opana::ICARUSBeamStructureAna::getRWMRelativeTime(int channel, double t) +{ + + if (fRWMTimes.empty()) + return icarus::timing::NoTime; + + auto rwm = fRWMTimes.at(channel); + if (!rwm.isValid()) + { + mf::LogTrace("ICARUSBeamStructureAna") << "No RWM signal for channel " << channel << " " + << "(Crate " << rwm.crate << ", Board " << rwm.digitizerLabel + << ", SpecialChannel " << rwm.specialChannel << ")" + << " in event " << m_event << " gate " << m_gate_name; + return icarus::timing::NoTime; + } + + double rwm_trigger = rwm.startTime; // rwm time w.r.t. trigger time [us] + return (t - rwm_trigger); +} + +// ----------------------------------------------------------------------------- + +double opana::ICARUSBeamStructureAna::getFlashBunchTime(std::vector channels, std::vector hit_rise_time_rwm) +{ + + double tfirst_left = std::numeric_limits::max(); + double tfirst_right = std::numeric_limits::max(); + + // if no RWM info available, all pmt_start_time_rwm are invalid + // return icarus::timing::NoTime as well for the flash + if (fRWMTimes.empty()) + return icarus::timing::NoTime; + + int nleft = 0; + int nright = 0; + for (std::size_t i = 0; i < hit_rise_time_rwm.size(); i++) + { + + int ch = channels[i]; + int side = getSideByChannel(ch); + double t = hit_rise_time_rwm[i]; // rise time w.r.t. rwm + + // if any RWM copy is missing (therefore missing for an entire PMT crate), + // it might not be possible to use the first hits (they might not have a RMW time) + // so return icarus::timing::NoTime as in other bad cases + if (!fRWMTimes[i].isValid()) + return icarus::timing::NoTime; + + // count hits separetely on the two walls + if (side == 0) + { + nleft++; + if (t < tfirst_left) + tfirst_left = t; + } + else if (side == 1) + { + nright++; + if (t < tfirst_right) + tfirst_right = t; + } + } + + // if there are no hits in one of the walls... very rare? + if (nleft < 1 || nright < 1) + { + mf::LogWarning("ICARUSBeamStructureAna") << "Flash " << m_flash_id << " doesn't have hits on both walls!" + << "Left: " << nleft << " t " << tfirst_left << " " + << "Right: " << nright << " t " << tfirst_right; + // return what we have... + return (tfirst_left < tfirst_right) ? tfirst_left : tfirst_right; + } + + return (tfirst_left + tfirst_right) / 2.; +} + +// ----------------------------------------------------------------------------- + +void opana::ICARUSBeamStructureAna::clear() +{ + + m_channel_id.clear(); + m_hit_start_time.clear(); + m_hit_peak_time.clear(); + m_hit_rise_time.clear(); + m_hit_start_time_rwm.clear(); + m_hit_peak_time_rwm.clear(); + m_hit_rise_time_rwm.clear(); + m_hit_pe.clear(); + + m_crthit_x.clear(); + m_crthit_y.clear(); + m_crthit_z.clear(); + m_crttime_us.clear(); + m_crtpmttimediff_ns.clear(); + m_crtsys.clear(); + m_crtregion.clear(); +} + +// ----------------------------------------------------------------------------- + +DEFINE_ART_MODULE(opana::ICARUSBeamStructureAna) diff --git a/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc b/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc index b47efb8f7..3943a5664 100644 --- a/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc +++ b/icaruscode/PMT/OpReco/ICARUSFlashAssAna_module.cc @@ -6,9 +6,10 @@ // Generated at Tue Jun 29 13:43:54 2021 by Andrea Scarpelli using cetskelgen // from cetlib version v3_11_01. // -// Module that dumps the association between Flashes and OpHit +// Module that dumps the association between Flashes and OpHit. +// These trees make up the optical information in the calibration ntuples. // -// mailto:ascarpel@bnl.gov +// mailto:ascarpel@bnl.gov, mvicenzi@bnl.gov //////////////////////////////////////////////////////////////////////// #include "art/Framework/Core/EDAnalyzer.h" @@ -40,217 +41,272 @@ #include "lardataobj/Simulation/BeamGateInfo.h" #include "lardataobj/RawData/TriggerData.h" #include "sbnobj/Common/Trigger/ExtraTriggerInfo.h" +#include "icaruscode/IcarusObj/PMTBeamSignal.h" +#include "sbnobj/ICARUS/PMT/Data/WaveformBaseline.h" #include "TTree.h" #include #include #include // std::accumulate +#include +#include - - -namespace opana { +namespace opana +{ class ICARUSFlashAssAna; } +class opana::ICARUSFlashAssAna : public art::EDAnalyzer +{ -class opana::ICARUSFlashAssAna : public art::EDAnalyzer { - - public: - - struct Config { +public: + struct Config + { - using Name = fhicl::Name; - using Comment = fhicl::Comment; + using Name = fhicl::Name; + using Comment = fhicl::Comment; - fhicl::Atom TriggerLabel { + fhicl::Atom TriggerLabel{ Name("TriggerLabel"), - Comment("Label for the Trigger fragment label") - }; + Comment("Label for the Trigger fragment label")}; - fhicl::Atom DumpWaveformsInfo { + fhicl::Atom DumpWaveformsInfo{ Name("DumpWaveformsInfo"), - Comment("Set the option to save some aggregated waveform information") - }; + Comment("Set the option to save some aggregated waveform information")}; - fhicl::Sequence OpDetWaveformLabels { + fhicl::Atom SaveRawWaveforms{ + Name("SaveRawWaveforms"), + Comment("Set to save the full raw::OpDetWaveforms")}; + + fhicl::Atom UseSharedBaseline{ + Name("UseSharedBaseline"), + Comment("Set the option to use icarus::WaveformBaseline")}; + + fhicl::Sequence OpDetWaveformLabels{ Name("OpDetWaveformLabels"), - Comment("Tags for the raw::OpDetWaveform data products") - }; + Comment("Tags for the raw::OpDetWaveform data products")}; - fhicl::Sequence OpHitLabels { + fhicl::Sequence BaselineLabels{ + Name("BaselineLabels"), + Comment("Tags for the icarus::WaveformBaseline data products")}; + + fhicl::Sequence OpHitLabels{ Name("OpHitLabels"), - Comment("Tags for the recob::OpHit data products") - }; + Comment("Tags for the recob::OpHit data products")}; - fhicl::Sequence FlashLabels { + fhicl::Sequence FlashLabels{ Name("FlashLabels"), - Comment("Tags for the recob::Flashe data products") - }; + Comment("Tags for the recob::Flashe data products")}; + + fhicl::Atom RWMLabel{ + Name("RWMLabel"), + Comment("Tag for the RWM std::vector data product")}; - fhicl::Atom PEOpHitThreshold { + fhicl::Atom PEOpHitThreshold{ Name("PEOpHitThreshold"), - Comment("Threshold in PE for an OpHit to be considered in the information calculated for a flash") - }; + Comment("Threshold in PE for an OpHit to be considered in the information calculated for a flash")}; - fhicl::Atom Debug { + fhicl::Atom Debug{ Name("Debug"), Comment("Be more verbose"), - false - }; - - }; - - using Parameters = art::EDAnalyzer::Table; - - explicit ICARUSFlashAssAna(Parameters const& config); - - ICARUSFlashAssAna(ICARUSFlashAssAna const&) = delete; - ICARUSFlashAssAna(ICARUSFlashAssAna&&) = delete; - ICARUSFlashAssAna& operator=(ICARUSFlashAssAna const&) = delete; - ICARUSFlashAssAna& operator=(ICARUSFlashAssAna&&) = delete; - - void analyze(art::Event const& e) override; - - void beginJob() override; - void endJob() override; - - template T Median( std::vector data ) const; - - geo::CryostatID::CryostatID_t getCryostatByChannel( int channel ); - - int getSideByChannel( const int channel ); - - void processOpHits( art::Event const& e, unsigned int cryo ); - - void processOpHitsFlash( std::vector> const &ophits, - int &multiplicity_left, int &multiplicity_right, - float &sum_pe_left, float &sum_pe_right, - float *xyz, - std::vector &pmt_start_time, - std::vector &pmt_pe, - std::vector &pmt_max_amplitude, - TTree *ophittree ); - - static std::string_view firstLine(std::string const& s, const char* endl = "\r"); - - private: - - art::InputTag fTriggerLabel; - bool fSaveWaveformInfo; - std::vector fOpDetWaveformLabels; - std::vector fOpHitLabels; - std::vector fFlashLabels; - float fPEOpHitThreshold; - bool fDebug; - - - TTree *fEventTree; - std::vector fOpDetWaveformTrees; - std::vector fOpFlashTrees; - std::vector fOpHitTrees; - std::vector fOpHitFlashTrees; - - int m_run; - int m_event; - int m_timestamp; - //int m_nflashes; - //int m_nophit; - short m_baseline; - short m_chargesum; - int m_nticks; - float m_beam_gate_start=-99999; - float m_beam_gate_width=-99999; - int m_beam_type=-1; - int m_trigger_type=-1; - unsigned int m_gate_type; - std::string m_gate_name; - uint64_t m_trigger_timestamp; - uint64_t m_gate_start_timestamp; - uint64_t m_trigger_gate_diff; - uint64_t lvdsCryoE[2]; - uint64_t lvdsCryoW[2]; - - int m_flash_id; - int m_multiplicity; - int m_multiplicity_left; - int m_multiplicity_right; - float m_sum_pe; - float m_sum_pe_left; - float m_sum_pe_right; - float m_flash_time; - //float m_flash_x; - //float m_flash_width_x; - float m_flash_y; - float m_flash_width_y; - float m_flash_z; - float m_flash_width_z; - std::vector m_pmt_time; - std::vector m_pmt_pe; - std::vector m_pmt_max_amplitude; - - int m_channel_id; - float m_integral; // in ADC x tick - float m_amplitude; // in ADC - float m_start_time; - float m_peak_time; - float m_rise_time; - float m_width; - float m_abs_start_time; - float m_pe; - float m_fast_to_total; - - std::vector m_pmt_x; - std::vector m_pmt_y; - std::vector m_pmt_z; - - geo::GeometryCore const* fGeom; - + false}; + }; + + using Parameters = art::EDAnalyzer::Table; + + explicit ICARUSFlashAssAna(Parameters const &config); + + ICARUSFlashAssAna(ICARUSFlashAssAna const &) = delete; + ICARUSFlashAssAna(ICARUSFlashAssAna &&) = delete; + ICARUSFlashAssAna &operator=(ICARUSFlashAssAna const &) = delete; + ICARUSFlashAssAna &operator=(ICARUSFlashAssAna &&) = delete; + + void analyze(art::Event const &e) override; + void beginJob() override; + + /// Compute median of waveform for baseline (fallback option) + template + T Median(std::vector data) const; + + /// Return cryostat from PMT channel_id + geo::CryostatID::CryostatID_t getCryostatByChannel(int channel); + + /// Return wall from PMT channel_id + int getSideByChannel(const int channel); + + /// Process OpHits in the absence of flashes + void processOpHits(art::Event const &e, unsigned int cryo); + + /// Process Ophits in the presence of flashes + void processOpHitsFlash(std::vector> const &ophits, + int &multiplicity_left, int &multiplicity_right, + float &sum_pe_left, float &sum_pe_right, + std::vector &pmt_start_time, + std::vector &pmt_rise_time, + std::vector &pmt_start_rwm_time, + std::vector &pmt_pe, + std::vector &pmt_amplitude, + TTree *ophittree); + + /// Return RWM-relative time from a trigger-relative time + float getRWMRelativeTime(int channel, float t); + + /// Return the RWM-relative flash interaction time + float getFlashBunchTime(std::vector pmt_start_time_rwm, + std::vector pmt_rise_time); + +private: + //---------- + // Input parameters + + art::InputTag fTriggerLabel; + bool fSaveWaveformInfo; + bool fSaveRawWaveforms; + bool fUseSharedBaseline; + std::vector fOpDetWaveformLabels; + std::vector fBaselineLabels; + std::vector fOpHitLabels; + std::vector fFlashLabels; + art::InputTag fRWMLabel; + float fPEOpHitThreshold; + bool fDebug; + + //---------- + // Output trees + + TTree *fEventTree; + std::vector fOpDetWaveformTrees; + std::vector fOpFlashTrees; + std::vector fOpHitTrees; + std::vector fOpHitFlashTrees; + + //---------------- + // Output variables + + // Common + int m_run; + int m_event; + int m_timestamp; + + // Event/trigger tree + float m_beam_gate_start = -99999; + float m_beam_gate_width = -99999; + int m_beam_type = -1; + int m_trigger_type = -1; + unsigned int m_gate_type; + std::string m_gate_name; + uint64_t m_trigger_timestamp; + uint64_t m_gate_start_timestamp; + uint64_t m_trigger_gate_diff; + uint64_t lvdsCryoE[2]; + uint64_t lvdsCryoW[2]; + uint16_t addersCryoE[2]; + uint16_t addersCryoW[2]; + + // Flash trees + int m_flash_id; + int m_multiplicity; + int m_multiplicity_left; + int m_multiplicity_right; + float m_sum_pe; + float m_sum_pe_left; + float m_sum_pe_right; + float m_flash_time; + float m_flash_time_rwm; + float m_flash_y; + float m_flash_width_y; + float m_flash_z; + float m_flash_width_z; + std::vector m_pmt_start_time; + std::vector m_pmt_rise_time; + std::vector m_pmt_start_time_rwm; + std::vector m_pmt_pe; + std::vector m_pmt_amplitude; + + // Ophit trees + int m_channel_id; + float m_integral; // in ADC x tick + float m_amplitude; // in ADC + float m_start_time; + float m_peak_time; + float m_rise_time; + float m_width; + float m_abs_start_time; + float m_start_time_rwm; + float m_peak_time_rwm; + float m_pe; + float m_fast_to_total; + + // Waveform trees + float m_wf_start; + short m_baseline; + short m_chargesum; + int m_nticks; + std::vector m_wf; + + // Geometry tree + std::vector m_pmt_x; + std::vector m_pmt_y; + std::vector m_pmt_z; + + //---------- + // Support variables/products + + geo::GeometryCore const *fGeom; + std::vector fRWMTimes; }; +// ---------------------------------------------------------------------------- + +opana::ICARUSFlashAssAna::ICARUSFlashAssAna(Parameters const &config) + : EDAnalyzer(config), + fTriggerLabel(config().TriggerLabel()), + fSaveWaveformInfo(config().DumpWaveformsInfo()), + fSaveRawWaveforms(config().SaveRawWaveforms()), + fUseSharedBaseline(config().UseSharedBaseline()), + fOpDetWaveformLabels(config().OpDetWaveformLabels()), + fBaselineLabels(config().BaselineLabels()), + fOpHitLabels(config().OpHitLabels()), + fFlashLabels(config().FlashLabels()), + fRWMLabel(config().RWMLabel()), + fPEOpHitThreshold(config().PEOpHitThreshold()), + fDebug(config().Debug()), + fGeom(lar::providerFrom()) +{ +} -opana::ICARUSFlashAssAna::ICARUSFlashAssAna(Parameters const& config) - : EDAnalyzer(config) - , fTriggerLabel( config().TriggerLabel() ) - , fSaveWaveformInfo( config().DumpWaveformsInfo() ) - , fOpDetWaveformLabels( config().OpDetWaveformLabels() ) - , fOpHitLabels( config().OpHitLabels() ) - , fFlashLabels( config().FlashLabels() ) - , fPEOpHitThreshold( config().PEOpHitThreshold() ) - , fDebug( config().Debug() ) - , fGeom( lar::providerFrom() ) -{ } - +// ---------------------------------------------------------------------------- -void opana::ICARUSFlashAssAna::beginJob() { +void opana::ICARUSFlashAssAna::beginJob() +{ art::ServiceHandle tfs; - TTree* fGeoTree = tfs->make("geotree", "geometry information" ); - fGeoTree->Branch("pmt_x",&m_pmt_x); - fGeoTree->Branch("pmt_y",&m_pmt_y); - fGeoTree->Branch("pmt_z",&m_pmt_z); - - for(size_t opch=0; opchNOpChannels(); ++opch) { + // Setting up the GEOMETRY tree + // Channel id corresponds to vector index + TTree *fGeoTree = tfs->make("geotree", "geometry information"); + fGeoTree->Branch("pmt_x", &m_pmt_x); + fGeoTree->Branch("pmt_y", &m_pmt_y); + fGeoTree->Branch("pmt_z", &m_pmt_z); - auto const PMTxyz = fGeom->OpDetGeoFromOpChannel(opch).GetCenter(); - - //std::cout << PMTxyz[0] << " " << PMTxyz[1] << " " << PMTxyz[2] << std::endl; + for (std::size_t opch = 0; opch < fGeom->NOpChannels(); ++opch) + { + auto const PMTxyz = fGeom->OpDetGeoFromOpChannel(opch).GetCenter(); m_pmt_x.push_back(PMTxyz.X()); m_pmt_y.push_back(PMTxyz.Y()); m_pmt_z.push_back(PMTxyz.Z()); - } fGeoTree->Fill(); - fEventTree = tfs->make("eventstree", "higher level information on the event" ); + // Settinp up the EVENT tree + // Trigger information and LVDS status + fEventTree = tfs->make("eventstree", "higher level information on the event"); fEventTree->Branch("run", &m_run, "run/I"); fEventTree->Branch("event", &m_event, "event/I"); fEventTree->Branch("timestamp", &m_timestamp, "timestamp/I"); - //fEventTree->Branch("nflashes", &m_nflashes, "nflashes/I"); - //fEventTree->Branch("nophits", &m_nophit, "nophits/I"); fEventTree->Branch("beam_gate_start", &m_beam_gate_start, "beam_gate_start/F"); fEventTree->Branch("beam_gate_width", &m_beam_gate_width, "beam_gate_width/F"); fEventTree->Branch("beam_type", &m_beam_type, "beam_type/I"); @@ -262,500 +318,634 @@ void opana::ICARUSFlashAssAna::beginJob() { fEventTree->Branch("trigger_gate_diff", &m_trigger_gate_diff, "trigger_gate_diff/l"); fEventTree->Branch("lvdsCryoE", &lvdsCryoE, "lvdsCryoE[2]/l"); fEventTree->Branch("lvdsCryoW", &lvdsCryoW, "lvdsCryoW[2]/l"); - + fEventTree->Branch("addersCryoE", &addersCryoE, "addersCryoE[2]/s"); + fEventTree->Branch("addersCryoW", &addersCryoW, "addersCryoW[2]/s"); + + // Setting up the WAVEFORM trees (one per product label) // This tree will hold some aggregated optical waveform information // The flag must be enabled to have the information saved - if( !fOpDetWaveformLabels.empty() && fSaveWaveformInfo ) { - - for( auto const & label : fOpDetWaveformLabels ) { + if (!fOpDetWaveformLabels.empty() && fSaveWaveformInfo) + { + + for (auto const &label : fOpDetWaveformLabels) + { - std::string name = label.label()+"wfttree"; + std::string name = label.label() + label.instance() + "wfttree"; std::string info = "TTree with aggregated optical waveform information with label: " + label.label(); - TTree* ttree = tfs->make(name.c_str(), info.c_str()); + TTree *ttree = tfs->make(name.c_str(), info.c_str()); ttree->Branch("run", &m_run, "run/I"); ttree->Branch("event", &m_event, "event/I"); ttree->Branch("timestamp", &m_timestamp, "timestamp/I"); ttree->Branch("channel_id", &m_channel_id, "channel_id/I"); + ttree->Branch("wf_start", &m_wf_start, "wf_start/F"); ttree->Branch("baseline", &m_baseline, "baseline/s"); ttree->Branch("chargesum", &m_chargesum, "chargesum/s"); ttree->Branch("nticks", &m_nticks, "nticks/I"); - - fOpDetWaveformTrees.push_back(ttree); + if (fSaveRawWaveforms) + ttree->Branch("wf", &m_wf); + fOpDetWaveformTrees.push_back(ttree); } - } - + // Setting up the OPHIT trees (one per cryostat) // This ttree will hold the ophit information when a flash is not found in the event // NB: information of the optical hits in events where flashes are present are lost - - for( auto const & label : fOpHitLabels ) { + for (auto const &label : fOpHitLabels) + { + + std::string name = label.label() + "_ttree"; + std::string info = "TTree for the recob::OpHit objects with label " + label.label() + " in events without flashes."; + + TTree *ttree = tfs->make(name.c_str(), info.c_str()); + ttree->Branch("run", &m_run, "run/I"); + ttree->Branch("event", &m_event, "event/I"); + ttree->Branch("timestamp", &m_timestamp, "timestamp/I"); + ttree->Branch("channel_id", &m_channel_id, "channel_id/I"); + ttree->Branch("integral", &m_integral, "integral/F"); + ttree->Branch("amplitude", &m_amplitude, "amplitude/F"); + ttree->Branch("start_time", &m_start_time, "start_time/F"); + ttree->Branch("peak_time", &m_peak_time, "peak_time/F"); + ttree->Branch("rise_time", &m_rise_time, "rise_time/F"); + ttree->Branch("abs_start_time", &m_abs_start_time, "abs_start_time/F"); + ttree->Branch("start_time_rwm", &m_start_time_rwm, "start_time_rwm/F"); + ttree->Branch("peak_time_rwm", &m_peak_time_rwm, "peak_time_rwm/F"); + ttree->Branch("pe", &m_pe, "pe/F"); + ttree->Branch("width", &m_width, "width/F"); + ttree->Branch("fast_to_total", &m_fast_to_total, "fast_to_total/F"); + + fOpHitTrees.push_back(ttree); + } + + // Setting up the OPFLASH/OPHITS trees (one per cryostat) + // These ttrees hold the information for the ophits and the flashes + // NB: information of the optical hits is stored differently when flashes are found + if (!fFlashLabels.empty()) + { - std::string name = label.label()+"_ttree"; - std::string info = "TTree for the recob::OpHit objects with label " + label.label() + " in events without flashes."; + for (auto const &label : fFlashLabels) + { - TTree* ttree = tfs->make(name.c_str(), info.c_str()); + // TTree for the flash in a given cryostat + std::string name = label.label() + "_flashtree"; + std::string info = "TTree for the recob::Flashes with label " + label.label(); + + TTree *ttree = tfs->make(name.c_str(), info.c_str()); ttree->Branch("run", &m_run, "run/I"); ttree->Branch("event", &m_event, "event/I"); ttree->Branch("timestamp", &m_timestamp, "timestamp/I"); - ttree->Branch("channel_id", &m_channel_id, "channel_id/I"); - ttree->Branch("integral", &m_integral, "integral/F"); - ttree->Branch("amplitude", &m_amplitude, "amplitude/F"); - ttree->Branch("start_time", &m_start_time, "start_time/F"); - ttree->Branch("peak_time", &m_peak_time, "peak_time/F"); - ttree->Branch("rise_time", &m_rise_time, "rise_time/F"); - ttree->Branch("abs_start_time", &m_abs_start_time, "abs_start_time/F"); - ttree->Branch("pe", &m_pe, "pe/F"); - ttree->Branch("width", &m_width, "width/F"); - ttree->Branch("fast_to_total", &m_fast_to_total, "fast_to_total/F"); - - fOpHitTrees.push_back(ttree); - - } - - - if ( !fFlashLabels.empty() ) { - - for( auto const & label : fFlashLabels ) { - - // TTree for the flash in a given cryostat - std::string name = label.label()+"_flashtree"; - std::string info = "TTree for the recob::Flashes with label "+label.label(); - - TTree* ttree = tfs->make(name.c_str(), info.c_str() ); - ttree->Branch("run", &m_run, "run/I"); - ttree->Branch("event", &m_event, "event/I"); - ttree->Branch("timestamp", &m_timestamp, "timestamp/I"); - ttree->Branch("flash_id", &m_flash_id, "flash_id/I"); - ttree->Branch("multiplicity", &m_multiplicity, "multiplicity/I"); - ttree->Branch("multiplicity_right", &m_multiplicity_right, "multiplicity_right/I" ); - ttree->Branch("multiplicity_left", &m_multiplicity_left, "multiplicity_left/I" ); - ttree->Branch("sum_pe", &m_sum_pe, "sum_pe/F"); - ttree->Branch("sum_pe_right", &m_sum_pe_right, "sum_pe_right/F"); - ttree->Branch("sum_pe_left", &m_sum_pe_left, "sum_pe_left/F"); - ttree->Branch("flash_time", &m_flash_time, "flash_time/F"); - //ttree->Branch("flash_x", &m_flash_x, "flash_x/F"); - //ttree->Branch("flash_width_x", &m_flash_width_x, "flash_width_x/F"); - ttree->Branch("flash_y", &m_flash_y, "flash_y/F"); - ttree->Branch("flash_width_y", &m_flash_width_y, "flash_width_y/F"); - ttree->Branch("flash_z", &m_flash_z, "flash_z/F"); - ttree->Branch("flash_width_z", &m_flash_width_z, "flash_width_z/F"); - ttree->Branch("pmt_x",&m_pmt_x); - ttree->Branch("pmt_y",&m_pmt_y); - ttree->Branch("pmt_z",&m_pmt_z); - ttree->Branch("time_pmt", & m_pmt_time); - ttree->Branch("pe_pmt", & m_pmt_pe ); - ttree->Branch("amplitude_pmt", &m_pmt_max_amplitude); - - fOpFlashTrees.push_back( ttree ); - - // Now the ttree for the OpHit associated in the flash - name = label.label()+"_ophittree"; - info = "Three for the recob::OpHit associated with an OpHitFlash"+label.label(); - - TTree* ophittree = tfs->make(name.c_str(), info.c_str() ); - ophittree->Branch("run", &m_run, "run/I"); - ophittree->Branch("event", &m_event, "event/I"); - ophittree->Branch("timestamp", &m_timestamp, "timestamp/I"); - ophittree->Branch("flash_id", &m_flash_id, "flash_id/I"); - ophittree->Branch("channel_id", &m_channel_id, "channel_id/I"); - ophittree->Branch("integral", &m_integral, "integral/F"); - ophittree->Branch("amplitude", &m_amplitude, "amplitude/F"); - ophittree->Branch("start_time", &m_start_time, "start_time/F"); - ophittree->Branch("peak_time", &m_peak_time, "peak_time/F"); - ophittree->Branch("rise_time", &m_rise_time, "rise_time/F"); - ophittree->Branch("abs_start_time", &m_abs_start_time, "abs_start_time/F"); - ophittree->Branch("pe", &m_pe, "pe/F"); - ophittree->Branch("width", &m_width, "width/F"); - ophittree->Branch("fast_to_total", &m_fast_to_total, "fast_to_total/F"); - - fOpHitFlashTrees.push_back( ophittree ); - + ttree->Branch("flash_id", &m_flash_id, "flash_id/I"); + ttree->Branch("multiplicity", &m_multiplicity, "multiplicity/I"); + ttree->Branch("multiplicity_right", &m_multiplicity_right, "multiplicity_right/I"); + ttree->Branch("multiplicity_left", &m_multiplicity_left, "multiplicity_left/I"); + ttree->Branch("sum_pe", &m_sum_pe, "sum_pe/F"); + ttree->Branch("sum_pe_right", &m_sum_pe_right, "sum_pe_right/F"); + ttree->Branch("sum_pe_left", &m_sum_pe_left, "sum_pe_left/F"); + ttree->Branch("flash_time", &m_flash_time, "flash_time/F"); + ttree->Branch("flash_time_rwm", &m_flash_time_rwm, "flash_time_rwm/F"); + ttree->Branch("flash_y", &m_flash_y, "flash_y/F"); + ttree->Branch("flash_width_y", &m_flash_width_y, "flash_width_y/F"); + ttree->Branch("flash_z", &m_flash_z, "flash_z/F"); + ttree->Branch("flash_width_z", &m_flash_width_z, "flash_width_z/F"); + ttree->Branch("pmt_x", &m_pmt_x); + ttree->Branch("pmt_y", &m_pmt_y); + ttree->Branch("pmt_z", &m_pmt_z); + ttree->Branch("time_pmt", &m_pmt_start_time); + ttree->Branch("time_pmt_rwm", &m_pmt_start_time_rwm); + ttree->Branch("pe_pmt", &m_pmt_pe); + ttree->Branch("amplitude_pmt", &m_pmt_amplitude); + + fOpFlashTrees.push_back(ttree); + + // Now the ttree for the OpHit associated in the flash + name = label.label() + "_ophittree"; + info = "Three for the recob::OpHit associated with an OpHitFlash" + label.label(); + + TTree *ophittree = tfs->make(name.c_str(), info.c_str()); + ophittree->Branch("run", &m_run, "run/I"); + ophittree->Branch("event", &m_event, "event/I"); + ophittree->Branch("timestamp", &m_timestamp, "timestamp/I"); + ophittree->Branch("flash_id", &m_flash_id, "flash_id/I"); + ophittree->Branch("channel_id", &m_channel_id, "channel_id/I"); + ophittree->Branch("integral", &m_integral, "integral/F"); + ophittree->Branch("amplitude", &m_amplitude, "amplitude/F"); + ophittree->Branch("start_time", &m_start_time, "start_time/F"); + ophittree->Branch("peak_time", &m_peak_time, "peak_time/F"); + ophittree->Branch("rise_time", &m_rise_time, "rise_time/F"); + ophittree->Branch("abs_start_time", &m_abs_start_time, "abs_start_time/F"); + ophittree->Branch("start_time_rwm", &m_start_time_rwm, "start_time_rwm/F"); + ophittree->Branch("peak_time_rwm", &m_peak_time_rwm, "peak_time_rwm/F"); + ophittree->Branch("pe", &m_pe, "pe/F"); + ophittree->Branch("width", &m_width, "width/F"); + ophittree->Branch("fast_to_total", &m_fast_to_total, "fast_to_total/F"); + + fOpHitFlashTrees.push_back(ophittree); } } - } +// ---------------------------------------------------------------------------- - -template - T opana::ICARUSFlashAssAna::Median( std::vector data ) const { - - std::nth_element( data.begin(), data.begin() + data.size()/2, data.end() ); - - return data[ data.size()/2 ]; +template +T opana::ICARUSFlashAssAna::Median(std::vector data) const +{ + std::nth_element(data.begin(), data.begin() + data.size() / 2, data.end()); + return data[data.size() / 2]; } +// ---------------------------------------------------------------------------- -geo::CryostatID::CryostatID_t opana::ICARUSFlashAssAna::getCryostatByChannel( int channel ) { - - - const geo::OpDetGeo& opdetgeo = fGeom->OpDetGeoFromOpChannel(channel); - geo::CryostatID::CryostatID_t cid = opdetgeo.ID().Cryostat ; +geo::CryostatID::CryostatID_t opana::ICARUSFlashAssAna::getCryostatByChannel(int channel) +{ + const geo::OpDetGeo &opdetgeo = fGeom->OpDetGeoFromOpChannel(channel); + geo::CryostatID::CryostatID_t cid = opdetgeo.ID().Cryostat; return cid; - } +// ---------------------------------------------------------------------------- -int opana::ICARUSFlashAssAna::getSideByChannel( const int channel ) { +int opana::ICARUSFlashAssAna::getSideByChannel(const int channel) +{ - /* + /* Channels are numbered from east to west, from North (cryo side) to South (beam side) - We look in the opposide direction wrt to the beam direction South->North: + We look in the opposide direction wrt to the beam direction South->North: - Left is the east wall of each cryostat; - - Right is the west side of each cryostat; - - - [ 0:89 ] and [180:269] are on the left, + - [ 0:89 ] and [180:269] are on the left, the return value of the function is 0; - - [ 90-179 ] and [ 270:359 ] are on the right, the return value of the function is 1; */ - int side = channel / 90; // always round down - return side % 2; } +// ---------------------------------------------------------------------------- -void opana::ICARUSFlashAssAna::processOpHits( art::Event const& e, unsigned int cryo ) { +float opana::ICARUSFlashAssAna::getRWMRelativeTime(int channel, float t) +{ + if (fRWMTimes.empty()) + return icarus::timing::NoTime; - if( fOpHitLabels.empty() ){ - - mf::LogError("ICARUSFlashAssAna") << "No recob::OpHit labels selected."; - - return; + auto rwm = fRWMTimes.at(channel); + if (!rwm.isValid()) + { + mf::LogTrace("ICARUSFlashAssAna") << "No RWM signal for channel " << channel << " " + << "(Crate " << rwm.crate << ", Board " << rwm.digitizerLabel + << ", SpecialChannel " << rwm.specialChannel << ")" + << " in event " << m_event << " gate " << m_gate_name; + return icarus::timing::NoTime; } - for( size_t iOpHitLabel=0; iOpHitLabel> ophit_handle; - e.getByLabel( label, ophit_handle ); - - - // We want our flashes to be valid and not empty - if( !ophit_handle.isValid() || ophit_handle->empty() ) { - mf::LogError("ICARUSFlashAssAna") - << "Invalid recob::OpHit with label '" << label.encode() << "'"; - continue; - } - - - for( auto const & ophit : *ophit_handle ) { + float rwm_trigger = rwm.startTime; // rwm time w.r.t. trigger time [us] + return (t - rwm_trigger); +} - //auto const & ophit = (*ophit_handle)[idx]; +// ---------------------------------------------------------------------------- + +float opana::ICARUSFlashAssAna::getFlashBunchTime(std::vector pmt_start_time_rwm, + std::vector pmt_rise_time) +{ + + float tfirst_left = std::numeric_limits::max(); + float tfirst_right = std::numeric_limits::max(); + + // if no RWM info available, all pmt_start_time_rwm are invalid + // return icarus::timing::NoTime as well for the flash + if (fRWMTimes.empty()) + return icarus::timing::NoTime; + + int nleft = 0; + int nright = 0; + for (std::size_t i = 0; i < pmt_start_time_rwm.size(); i++) + { + + int side = getSideByChannel(i); + float t = pmt_start_time_rwm[i] + pmt_rise_time[i]; // rise time w.r.t. rwm + + // exclude channels that have no signals + // these have time == 0 + if (t > -1e-10 && t < 1e-10) + continue; + + // if any RWM copy is missing (therefore missing for an entire PMT crate), + // it might not be possible to use the first hits (they might not have a RMW time) + // so return icarus::timing::NoTime as in other bad cases + if (!fRWMTimes[i].isValid()) + return icarus::timing::NoTime; + + // count hits separetely on the two walls + if (side == 0) + { + nleft++; + if (t < tfirst_left) + tfirst_left = t; + } + else if (side == 1) + { + nright++; + if (t < tfirst_right) + tfirst_right = t; + } + } - const int channel_id = ophit.OpChannel(); + // if there are no hits in one of the walls... + if (nleft < 1 || nright < 1) + { + mf::LogWarning("ICARUSFlashAssAna") << "Flash " << m_flash_id << " doesn't have hits on both walls!" + << "Left: " << nleft << " t " << tfirst_left << " " + << "Right: " << nright << " t " << tfirst_right; + // return what we have... + return (tfirst_left < tfirst_right) ? tfirst_left : tfirst_right; + } - if( getCryostatByChannel(channel_id) != cryo ){ continue; } + return (tfirst_left + tfirst_right) / 2.; +} - m_channel_id = channel_id; - m_integral = ophit.Area(); // in ADC x tick - m_amplitude = ophit.Amplitude(); // in ADC - m_start_time = ophit.StartTime(); - m_peak_time = ophit.PeakTime(); - m_rise_time = ophit.RiseTime(); - m_width = ophit.Width(); - m_abs_start_time = ophit.PeakTimeAbs() + (m_start_time - m_peak_time); - m_pe = ophit.PE(); - m_fast_to_total = ophit.FastToTotal(); +// ---------------------------------------------------------------------------- - fOpHitTrees[iOpHitLabel]->Fill(); +void opana::ICARUSFlashAssAna::processOpHits(art::Event const &e, unsigned int cryo) +{ - } + // if no OpHits have been selected at all (and no flashes as well!) + if (fOpHitLabels.empty()) + { + mf::LogError("ICARUSFlashAssAna") << "No recob::OpHit labels selected."; + return; } - return; - -} + for (std::size_t iOpHitLabel = 0; iOpHitLabel < fOpHitLabels.size(); iOpHitLabel++) + { + auto const label = fOpHitLabels[iOpHitLabel]; + auto const &ophits = e.getProduct>(label); -void opana::ICARUSFlashAssAna::processOpHitsFlash( std::vector> const &ophits, - int &multiplicity_left, int &multiplicity_right, - float &sum_pe_left, float &sum_pe_right, - float *xyz, - std::vector &pmt_start_time, - std::vector &pmt_pe, - std::vector &pmt_max_amplitude, - TTree *ophittree ) { - + // we want our ophits to be valid and not empty + if (ophits.empty()) + { + mf::LogError("ICARUSFlashAssAna") << "Invalid recob::OpHit with label '" << label.encode() << "'"; + continue; + } - std::unordered_map sumpe_map; + for (auto const &ophit : ophits) + { - // We caluclate the total charge clustered in the flash per channel taking part to the flash - for( auto const ophit : ophits ) { + const int channel_id = ophit.OpChannel(); + if (getCryostatByChannel(channel_id) != cryo) + { + continue; + } - if ( ophit->PE() < fPEOpHitThreshold ) { continue; } + m_channel_id = channel_id; + m_integral = ophit.Area(); // in ADC x tick + m_amplitude = ophit.Amplitude(); // in ADC + m_width = ophit.Width(); + m_pe = ophit.PE(); + m_fast_to_total = ophit.FastToTotal(); + + // save times: start, peak, rise + // rise is relative to start + m_start_time = ophit.StartTime(); + m_peak_time = ophit.PeakTime(); + m_rise_time = ophit.RiseTime(); + m_start_time_rwm = getRWMRelativeTime(channel_id, m_start_time); + m_peak_time_rwm = getRWMRelativeTime(channel_id, m_peak_time); + m_abs_start_time = ophit.PeakTimeAbs() + (m_start_time - m_peak_time); + + fOpHitTrees[iOpHitLabel]->Fill(); + } + } +} - const int channel_id = ophit->OpChannel(); +// ---------------------------------------------------------------------------- + +void opana::ICARUSFlashAssAna::processOpHitsFlash(std::vector> const &ophits, + int &multiplicity_left, int &multiplicity_right, + float &sum_pe_left, float &sum_pe_right, + std::vector &pmt_start_time, + std::vector &pmt_rise_time, + std::vector &pmt_start_time_rwm, + std::vector &pmt_pe, + std::vector &pmt_amplitude, + TTree *ophittree) +{ + + // we calculate the total charge clustered in the flash per channel taking part to the flash + // at the same time we store the times of the first ophit in each channel and its amplitude + // same loop is used to saved OPHITS info + + std::unordered_map sumpe_map; + for (auto const ophit : ophits) + { + + if (ophit->PE() < fPEOpHitThreshold) + { + continue; + } - sumpe_map[ channel_id ]+=ophit->PE() ; + const int channel_id = ophit->OpChannel(); + sumpe_map[channel_id] += ophit->PE(); m_channel_id = channel_id; - m_integral = ophit->Area(); // in ADC x tick + m_integral = ophit->Area(); // in ADC x tick m_amplitude = ophit->Amplitude(); // in ADC + m_width = ophit->Width(); + m_pe = ophit->PE(); + m_fast_to_total = ophit->FastToTotal(); + + // save times: start, peak, rise + // rise is relative to start m_start_time = ophit->StartTime(); m_peak_time = ophit->PeakTime(); m_rise_time = ophit->RiseTime(); - m_width = ophit->Width(); + m_start_time_rwm = getRWMRelativeTime(channel_id, m_start_time); + m_peak_time_rwm = getRWMRelativeTime(channel_id, m_peak_time); m_abs_start_time = ophit->PeakTimeAbs() + (m_start_time - m_peak_time); - m_pe = ophit->PE(); - m_fast_to_total = ophit->FastToTotal(); pmt_pe[channel_id] += ophit->PE(); - if( ( pmt_start_time[channel_id] == 0 ) || ( pmt_start_time[channel_id] > m_start_time )) { + // select the first ophit (by time) in each channel + if ((pmt_start_time[channel_id] == 0) || (pmt_start_time[channel_id] > m_start_time)) + { pmt_start_time[channel_id] = m_start_time; - pmt_max_amplitude[channel_id] = m_amplitude; + pmt_rise_time[channel_id] = m_rise_time; + pmt_start_time_rwm[channel_id] = m_start_time_rwm; + pmt_amplitude[channel_id] = m_amplitude; } - ophittree->Fill(); - } - m_multiplicity_left = std::accumulate( sumpe_map.begin(), sumpe_map.end(), 0, - [&](int value, const std::map::value_type& p) { - return getSideByChannel(p.first)==0 ? ++value : value ; - }); - - m_multiplicity_right =std::accumulate( sumpe_map.begin(), sumpe_map.end(), 0, - [&](int value, const std::map::value_type& p) { - return getSideByChannel(p.first)==1 ? ++value : value ; - }); - - m_sum_pe_left = std::accumulate( sumpe_map.begin(), sumpe_map.end(), 0.0, - [&](float value, const std::map::value_type& p) { - return getSideByChannel(p.first)==0 ? value+p.second : value ; - }); - - m_sum_pe_right = std::accumulate( sumpe_map.begin(), sumpe_map.end(), 0.0, - [&](float value, const std::map::value_type& p) { - return getSideByChannel(p.first)==1 ? value+p.second : value ; - }); - - //for( int i=0; i<3; i++ ){ xyz[i] /= (m_sum_pe_left+ m_sum_pe_right); } - -} + m_multiplicity_left = std::accumulate(sumpe_map.begin(), sumpe_map.end(), 0, + [&](int value, const std::map::value_type &p) + { return getSideByChannel(p.first) == 0 ? ++value : value; }); + m_multiplicity_right = std::accumulate(sumpe_map.begin(), sumpe_map.end(), 0, + [&](int value, const std::map::value_type &p) + { return getSideByChannel(p.first) == 1 ? ++value : value; }); -void opana::ICARUSFlashAssAna::endJob() { + m_sum_pe_left = std::accumulate(sumpe_map.begin(), sumpe_map.end(), 0.0, + [&](float value, const std::map::value_type &p) + { return getSideByChannel(p.first) == 0 ? value + p.second : value; }); + m_sum_pe_right = std::accumulate(sumpe_map.begin(), sumpe_map.end(), 0.0, + [&](float value, const std::map::value_type &p) + { return getSideByChannel(p.first) == 1 ? value + p.second : value; }); } +// ---------------------------------------------------------------------------- -void opana::ICARUSFlashAssAna::analyze(art::Event const& e) { - +void opana::ICARUSFlashAssAna::analyze(art::Event const &e) +{ + // Collect global event metadata m_run = e.id().run(); m_event = e.id().event(); - m_timestamp = e.time().timeHigh(); // precision to the second - - - /* - This part is for the trigger information - */ + m_timestamp = e.time().timeHigh(); // precision to the second + // ----- + // TRIGGER INFO + // We work out the trigger information here - // We work out the trigger information here - if( !fTriggerLabel.empty() ) { + if (!fTriggerLabel.empty()) + { - // Beam information - art::Handle> beamgate_handle; - e.getByLabel( fTriggerLabel, beamgate_handle ); - - if( beamgate_handle.isValid() ) { + // Beam gate information + auto const &beamgateInfo = e.getProduct>(fTriggerLabel); - for( auto const & beamgate : *beamgate_handle ) { + if (!beamgateInfo.empty()) + { + for (auto const &beamgate : beamgateInfo) + { + m_beam_gate_start = beamgate.Start(); + m_beam_gate_width = beamgate.Width(); + m_beam_type = beamgate.BeamType(); + } + } + else + { + mf::LogError("ICARUSFlashAssAna") << "No sim::BeamGateInfo associated to label: " << fTriggerLabel.label() << "\n"; + } - m_beam_gate_start = beamgate.Start(); - m_beam_gate_width = beamgate.Width(); - m_beam_type = beamgate.BeamType() ; + // Now trigger information + auto const &extraInfo = e.getProduct(fTriggerLabel); - } + if (extraInfo.isValid()) + { - } - - else { - mf::LogError("ICARUSFlashAssAna") << "No sim::BeamGateInfo associated to label: " << fTriggerLabel.label() << "\n" ; - } + sbn::triggerSource bit = extraInfo.sourceType; + m_gate_type = (unsigned int)bit; + m_gate_name = bitName(bit); - // Now trigger information - art::Handle trigger_handle; - e.getByLabel( fTriggerLabel, trigger_handle ); + m_trigger_type = value(extraInfo.triggerType); + m_trigger_timestamp = extraInfo.triggerTimestamp; + m_gate_start_timestamp = extraInfo.beamGateTimestamp; + m_trigger_gate_diff = m_trigger_timestamp - m_gate_start_timestamp; - if( trigger_handle.isValid() ) { + // majority lvds info + lvdsCryoE[0] = extraInfo.cryostats[0].LVDSstatus[0]; + lvdsCryoE[1] = extraInfo.cryostats[0].LVDSstatus[1]; + lvdsCryoW[0] = extraInfo.cryostats[1].LVDSstatus[0]; + lvdsCryoW[1] = extraInfo.cryostats[1].LVDSstatus[1]; - sbn::triggerSource bit = trigger_handle->sourceType; + // adders lvds info + addersCryoE[0] = extraInfo.cryostats[0].sectorStatus[0]; + addersCryoE[1] = extraInfo.cryostats[0].sectorStatus[1]; + addersCryoW[0] = extraInfo.cryostats[1].sectorStatus[0]; + addersCryoW[1] = extraInfo.cryostats[1].sectorStatus[1]; + } + else + { + mf::LogError("ICARUSFlashAssAna") << "No raw::Trigger associated to label: " << fTriggerLabel.label() << "!"; + } + } - m_gate_type = (unsigned int)bit; - m_gate_name = bitName(bit); - m_trigger_type = value( trigger_handle->triggerType ); - m_trigger_timestamp = trigger_handle->triggerTimestamp; - m_gate_start_timestamp = trigger_handle->beamGateTimestamp; - m_trigger_gate_diff = trigger_handle->triggerTimestamp - trigger_handle->beamGateTimestamp; - lvdsCryoE[0] = trigger_handle->cryostats[0].LVDSstatus[0]; - lvdsCryoE[1] = trigger_handle->cryostats[0].LVDSstatus[1]; - lvdsCryoW[0] = trigger_handle->cryostats[1].LVDSstatus[0]; - lvdsCryoW[1] = trigger_handle->cryostats[1].LVDSstatus[1]; + // ----- + // RWM INFO + // We work out the RWM information here + // it might be empty if offbeam or missing, bu that's okay! - } - else{ - mf::LogError("ICARUSFlashAssAna") << "No raw::Trigger associated to label: " << fTriggerLabel.label() << "\n" ; - } + if (!fRWMLabel.empty()) + { + fRWMTimes = e.getProduct>(fRWMLabel); + if (fRWMTimes.empty()) + mf::LogTrace("ICARUSFlashAssAna") << "Data product std::vector for '" + << fRWMLabel.label() << "' is empty in " << m_gate_name << " event!"; } - else { - mf::LogError("ICARUSFlashAssAna") << "Trigger Data product " << fTriggerLabel.label() << " not found!\n" ; - } + // ----- + // WAVEFORM INFO + // Now we work on the waveforms if we are allowed to + // Full raw waveforms are dumped only if option is set explicitly! + // Waveforms are a complex business: if running at stage0, they're all available + // if running at stage1, only on-beam waveforms are available but it's a different product! + if (!fOpDetWaveformLabels.empty() && fSaveWaveformInfo) + { - /* - Now we work on the waveforms if we are allowed to - */ + for (std::size_t i = 0; i < fOpDetWaveformLabels.size(); i++) + { - if( !fOpDetWaveformLabels.empty() && fSaveWaveformInfo ) { - - for( size_t i=0; i>(wflabel); + auto const &baselines = e.getProduct>(bslabel); - auto const label = fOpDetWaveformLabels[i]; + // waveforms data product ("daqPMT") is dropped after stage0, and only "daqPMTonbeam" is kept + // check if collection is present and print warning otherwise (requires change in fhicl config) + if (!waveforms.empty()) + { - art::Handle> wfm_handle; - e.getByLabel( label, wfm_handle ); + std::size_t idx = 0; + for (auto const &wave : waveforms) + { - if( wfm_handle.isValid() && !wfm_handle->empty() ) { + m_channel_id = wave.ChannelNumber(); + m_nticks = wave.Waveform().size(); + m_wf_start = wave.TimeStamp(); - for( auto const & wave : *wfm_handle ){ + // try using icarus::WaveformBaseline, default to median if not + // doesn't work with daqPMTonbeam as order gets messed up, so check if size is the same + // FIXME: move to use art:Assn which can work both for daqPMT and daqPMTonbeam? + m_baseline = (fUseSharedBaseline && (baselines.size() == waveforms.size())) ? baselines.at(idx).baseline() : Median(wave.Waveform()); - m_channel_id = wave.ChannelNumber(); - m_nticks = wave.Waveform().size(); - m_baseline = Median( wave.Waveform() ); - m_chargesum = std::accumulate( wave.Waveform().begin(), wave.Waveform().end(), 0, - [ & ](short x, short y){ return ((m_baseline-x) + (m_baseline-y)) ; } ); + m_chargesum = std::accumulate(wave.Waveform().begin(), wave.Waveform().end(), 0, + [&](short x, short y) + { return ((m_baseline - x) + (m_baseline - y)); }); - fOpDetWaveformTrees[i]->Fill(); - - } + // if required, save the full waveforms as well + if (fSaveRawWaveforms) + m_wf = wave.Waveform(); + + fOpDetWaveformTrees[i]->Fill(); + idx++; } } + else + { + mf::LogWarning("ICARUSFlashAssAna") << "Data product std::vector for " << wflabel.label() + << " is missing! Was it dropped before this stage?"; + } } + } + // ----- + // FLASHES/OPHITS INFO + // Now we take care of the flashes: + // we separate the case where we have a flash and the case where we don't have a flash + // the difference is in how ophits are stored... - /* - Now we take care of the flashes: we separate the case where we have a flash and the case where we have not a flash - */ - - if ( !fFlashLabels.empty() ) { + if (!fFlashLabels.empty()) + { - // Hold the cryostat information - std::vector cids; - + // hold the cryostat information + std::vector cids; - for ( size_t iFlashLabel=0; iFlashLabel>(label); + auto const &flashes = *flash_handle; - art::Handle> flash_handle; - e.getByLabel( label, flash_handle ); - - // We want our flashes to be valid and not empty - if( !flash_handle.isValid() ) { - mf::LogError("ICARUSFlashAssAna") - << "Not found a recob::OpFlash with label '" << label.encode() << "'"; - } else if ( flash_handle->empty() ) { + // want our flashes to be valid and not empty + if (flashes.empty()) + { mf::LogWarning("ICARUSFlashAssAna") - << "No recob::OpFlash in collection with label '" << label.encode() << "'"; + << "No recob::OpFlash in collection with label '" << label.encode() << "'"; } - else { + else + { - art::FindManyP ophitsPtr( flash_handle, e, label ); + art::FindManyP ophitsPtr(flash_handle, e, label); - for ( size_t idx=0; idxsize(); idx++ ) { + std::size_t idx = 0; + for (auto const &flash : flashes) + { m_pmt_pe.resize(360); - m_pmt_time.resize(360); - m_pmt_max_amplitude.resize(360); + m_pmt_start_time.resize(360); + m_pmt_rise_time.resize(360); + m_pmt_start_time_rwm.resize(360); + m_pmt_amplitude.resize(360); m_flash_id = idx; - auto const & flash = (*flash_handle)[idx]; - m_flash_time = flash.Time(); m_sum_pe = flash.TotalPE(); - - auto const & ophits = ophitsPtr.at(idx); + m_flash_y = flash.YCenter(); + m_flash_width_y = flash.YWidth(); + m_flash_z = flash.ZCenter(); + m_flash_width_z = flash.ZWidth(); - // We keep track of the cryistats where the flashes are found; + auto const &ophits = ophitsPtr.at(idx); + + // we keep track of the cryistats where the flashes are found; geo::CryostatID::CryostatID_t cid = getCryostatByChannel(ophits.front()->OpChannel()); - auto const found = std::find(cids.begin(), cids.end(), cid); - if( found != cids.end() ){ - cids.push_back( cid ); + if (found != cids.end()) + { + cids.push_back(cid); } - // Get the multiplicity, the position and the number of PE per Side - float xyz[3] = {0.0, 0.0, 0.0}; - processOpHitsFlash( ophits, - m_multiplicity_left, m_multiplicity_right, - m_sum_pe_left, m_sum_pe_right, xyz, - m_pmt_time, m_pmt_pe, m_pmt_max_amplitude, - fOpHitFlashTrees[iFlashLabel] ); + // get the multiplicity, the number of PE per side, and the RWM-relative flash time + // also store the first ophits for every channel in the flash + processOpHitsFlash(ophits, + m_multiplicity_left, m_multiplicity_right, + m_sum_pe_left, m_sum_pe_right, m_pmt_start_time, + m_pmt_rise_time, m_pmt_start_time_rwm, m_pmt_pe, m_pmt_amplitude, + fOpHitFlashTrees[iFlashLabel]); - m_multiplicity = m_multiplicity_left+m_multiplicity_right; + m_multiplicity = m_multiplicity_left + m_multiplicity_right; - //m_flash_x = 0.0; - //m_flash_width_x = 0.0; - m_flash_y = flash.YCenter(); - m_flash_width_y = flash.YWidth(); - m_flash_z = flash.ZCenter(); - m_flash_width_z = flash.ZWidth(); + // get the flash interaction time w.r.t. RWM + // this is currently the mean between the first ophits on opposite walls + m_flash_time_rwm = getFlashBunchTime(m_pmt_start_time_rwm, m_pmt_rise_time); fOpFlashTrees[iFlashLabel]->Fill(); m_pmt_pe.clear(); - m_pmt_time.clear(); + m_pmt_start_time.clear(); + m_pmt_rise_time.clear(); + m_pmt_start_time_rwm.clear(); + idx++; } } - } + } - // If the flashes did not cover all three cryostats.. + // If the flashes did not cover all three cryostats.. // ..well, we save the ophits on what is missing - for( unsigned int cid=0; cidNcryostats(); cid++ ){ + for (unsigned int cid = 0; cid < fGeom->Ncryostats(); cid++) + { - auto const found = std::find( cids.begin(), cids.end(), cid ); - if( found == cids.end() ){ - processOpHits(e, cid); - } + auto const found = std::find(cids.begin(), cids.end(), cid); + if (found == cids.end()) + { + processOpHits(e, cid); + } } } - else { - - mf::LogError("ICARUSFlashAssAna") - << "No recob::OpFlash labels selected\n"; + else + { + mf::LogError("ICARUSFlashAssAna") << "No recob::OpFlash labels selected"; - // We save the ophits anyways even in absence of flashes - for( unsigned int cid=0; cidNcryostats(); cid++ ){ - processOpHits(e, cid); + // we save the ophits anyways even in absence of flashes + for (unsigned int cid = 0; cid < fGeom->Ncryostats(); cid++) + { + processOpHits(e, cid); } - } - fEventTree->Fill(); - } - DEFINE_ART_MODULE(opana::ICARUSFlashAssAna) diff --git a/icaruscode/PMT/OpReco/driver/run_beamstructure_ana.fcl b/icaruscode/PMT/OpReco/driver/run_beamstructure_ana.fcl new file mode 100644 index 000000000..f279723b9 --- /dev/null +++ b/icaruscode/PMT/OpReco/driver/run_beamstructure_ana.fcl @@ -0,0 +1,82 @@ +################################################################# +# File: run_beamstructure_ana.fcl +# Author: M. Vicenzi (mvicenzi@bnl.gov) +# +# Description: +# Simple off-the-shelf example to run the beam structure analysis. +# It produces all ingredients for the light-only reconstruction of the beam bunch structure. +# This is meant to be run directly on raw files as it defines all the needed producers. +# Output contains both RWM/EW information and OpFlash/CRT timing information. + +#include "services_common_icarus.fcl" +#include "channelmapping_icarus.fcl" +#include "timing_icarus.fcl" + +#include "decoderdefs_icarus.fcl" +#include "icarus_ophitfinder.fcl" +#include "icarus_flashfinder.fcl" +#include "icarus_opana_modules.fcl" +#include "timing_beam.fcl" + +#include "crt_decoderdefs_icarus.fcl" +#include "crthitproducer.fcl" +#include "crtpmtmatchingproducer.fcl" +#include "icarus_FilterCRTPMTMatching.fcl" + +process_name: beamana + +services: +{ + @table::icarus_art_services # from services_common_icarus.fcl + @table::icarus_basic_services # from services_basic_icarus.fcl + @table::icarus_geometry_services + DetectorClocksService: @local::icarus_detectorclocks + IICARUSChannelMap: @local::icarus_channelmappinggservice + IPMTTimingCorrectionService: @local::icarus_pmttimingservice + TFileService: { fileName: "supplemental-%ifb-%p.root" } +} + +physics: +{ + producers: + { + triggerconfig: @local::extractTriggerConfig + pmtconfig: @local::extractPMTconfig + + daqTrigger: @local::decodeTriggerAutodetect + daqPMT: @local::decodePMT + daqCRT: @local::crtdaq_icarus + + ophituncorrected: @local::icarus_ophit_data + ophit: @local::icarus_ophit_timing_correction + beamTiming: @local::icarus_beam_signal_extractor + opflashCryoE: @local::ICARUSSimpleFlashDataCryoE + opflashCryoW: @local::ICARUSSimpleFlashDataCryoW + + crthit: @local::standard_crthitproducer + crtpmt: @local::standard_crtpmtmatchingproducer + } + + analyzers: + { + beamana: @local::ICARUSBeamStructureAna + } + + my_producers: [ triggerconfig, pmtconfig, daqTrigger, daqPMT, daqCRT, ophituncorrected, ophit, beamTiming, opflashCryoE, opflashCryoW, crthit, crtpmt ] + my_analyzers: [ beamana ] + + trigger_paths: [ my_producers ] + end_paths: [ my_analyzers ] +} + +### REQUIRED PRODUCERS ### +physics.producers.daqTrigger.DecoderTool.Decoders[0].ToolConfig.TrigConfigLabel: triggerconfig +physics.producers.daqTrigger.DecoderTool.Decoders[1].ToolConfig.TrigConfigLabel: triggerconfig +physics.producers.daqTrigger.DecoderTool.Decoders[2].ToolConfig.TrigConfigLabel: triggerconfig + +physics.producers.daqPMT.PMTconfigTag: pmtconfig # required +physics.producers.daqPMT.TriggerTag: daqTrigger # required +physics.producers.ophit.InputLabels: [ "ophituncorrected" ] + +physics.producers.beamTiming.DebugTrees: true +physics.analyzers.beamana.RWMLabel: "beamTiming:RWM" diff --git a/icaruscode/PMT/OpReco/driver/run_flashana.fcl b/icaruscode/PMT/OpReco/driver/run_flashana.fcl new file mode 100644 index 000000000..0e26d3d4f --- /dev/null +++ b/icaruscode/PMT/OpReco/driver/run_flashana.fcl @@ -0,0 +1,74 @@ +####################################################### +# File: run_flashana.fcl +# Author: Matteo Vicenzi (mvicenzi@bnl.gov) +# Description: +# Simple off-the-shelf example to run the light analysis +# directly from raw files + +#include "services_common_icarus.fcl" +#include "channelmapping_icarus.fcl" +#include "timing_icarus.fcl" +#include "timing_beam.fcl" + +#include "decoderdefs_icarus.fcl" +#include "icarus_ophitfinder.fcl" +#include "icarus_flashfinder.fcl" +#include "icarus_opana_modules.fcl" + +process_name: flashana + +services: +{ + @table::icarus_art_services # from services_common_icarus.fcl + @table::icarus_basic_services # from services_basic_icarus.fcl + @table::icarus_geometry_services + DetectorClocksService: @local::icarus_detectorclocks + IICARUSChannelMap: @local::icarus_channelmappinggservice + IPMTTimingCorrectionService: @local::icarus_pmttimingservice + TFileService: { fileName: "supplemental-%ifb-%p.root" } +} + +physics: +{ + producers: + { + triggerconfig: @local::extractTriggerConfig + pmtconfig: @local::extractPMTconfig + + daqTrigger: @local::decodeTriggerAutodetect + daqPMT: @local::decodePMT + + pmtbaselines: @local::icarus_opreco_pedestal_fromchannel_data # from icarus_ophitfinder.fcl + ophituncorrected: @local::icarus_ophit_data + ophit: @local::icarus_ophit_timing_correction + beamTiming: @local::icarus_beam_signal_extractor + opflashCryoE: @local::ICARUSSimpleFlashDataCryoE + opflashCryoW: @local::ICARUSSimpleFlashDataCryoW + } + + analyzers: + { + flashana: @local::ICARUSFlashAssAna + } + + my_producers: [ triggerconfig, pmtconfig, daqTrigger, daqPMT, pmtbaselines, ophituncorrected, ophit, beamTiming, opflashCryoE, opflashCryoW ] + my_analyzers: [ flashana ] + + trigger_paths: [ my_producers ] + end_paths: [ my_analyzers ] +} + +### REQUIRED PRODUCERS ### +physics.producers.daqTrigger.DecoderTool.Decoders[0].ToolConfig.TrigConfigLabel: triggerconfig +physics.producers.daqTrigger.DecoderTool.Decoders[1].ToolConfig.TrigConfigLabel: triggerconfig +physics.producers.daqTrigger.DecoderTool.Decoders[2].ToolConfig.TrigConfigLabel: triggerconfig + +physics.producers.daqPMT.PMTconfigTag: pmtconfig # required +physics.producers.daqPMT.TriggerTag: daqTrigger # required + +physics.producers.ophit.InputLabels: [ "ophituncorrected" ] +physics.producers.beamTiming.DebugTrees: false + +physics.analyzers.flashana.DumpWaveformsInfo: true +physics.analyzers.flashana.SaveRawWaveforms: false +physics.analyzers.flashana.UseSharedBaseline: true diff --git a/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl b/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl index 255b82170..e64474a00 100644 --- a/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl +++ b/icaruscode/PMT/OpReco/fcl/icarus_opana_modules.fcl @@ -104,11 +104,24 @@ ICARUSFlashAssAna: { module_type: "ICARUSFlashAssAna" TriggerLabel: "daqTrigger" DumpWaveformsInfo: true + SaveRawWaveforms: false + UseSharedBaseline: true OpDetWaveformLabels: ["daqPMT"] + BaselineLabels: ["pmtbaselines"] OpHitLabels: ["ophit"] FlashLabels: ["opflashCryoE", "opflashCryoW"] + RWMLabel: "beamTiming:RWM" PEOpHitThreshold: 0 Debug: false } +ICARUSBeamStructureAna: { + module_type: "ICARUSBeamStructureAna" + FlashLabels: ["opflashCryoE", "opflashCryoW"] + TriggerLabel: "daqTrigger" + RWMLabel: "beamTiming:RWM" + TriggerConfigLabel: "triggerconfig" + CRTPMTMatchingLabel: "crtpmt" +} + END_PROLOG diff --git a/icaruscode/Timing/CMakeLists.txt b/icaruscode/Timing/CMakeLists.txt index 0e652408e..85fc1cc05 100644 --- a/icaruscode/Timing/CMakeLists.txt +++ b/icaruscode/Timing/CMakeLists.txt @@ -2,8 +2,13 @@ cet_enable_asserts() set( MODULE_LIBRARIES lardataobj::RecoBase + sbnobj::Common_Trigger larcore::Geometry_Geometry_service - lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + lardata::DetectorClocksService + art_root_io::TFileService_service + art_root_io::tfile_support + lardataobj::RawData + ROOT::Tree ) set( LIB_LIBRARIES art::Framework_Services_Registry @@ -19,7 +24,7 @@ set( SERVICE_LIBRARIES icaruscode_Timing icaruscode_IcarusObj larcore::Geometry_Geometry_service - lardata::DetectorInfoServices_DetectorClocksServiceStandard_service + lardata::DetectorClocksService ) file(GLOB lib_srcs *.cxx) @@ -30,6 +35,8 @@ cet_build_plugin( OpHitTimingCorrection art::module LIBRARIES PUBLIC ${MODULE_LI cet_build_plugin( PMTTimingCorrectionService art::service LIBRARIES PUBLIC ${SERVICE_LIBRARIES}) +cet_build_plugin( PMTBeamSignalsExtractor art::producer LIBRARIES PUBLIC ${MODULE_LIBRARIES}) + install_headers() install_fhicl() install_source() diff --git a/icaruscode/Timing/PMTBeamSignalsExtractor_module.cc b/icaruscode/Timing/PMTBeamSignalsExtractor_module.cc new file mode 100644 index 000000000..c190f770f --- /dev/null +++ b/icaruscode/Timing/PMTBeamSignalsExtractor_module.cc @@ -0,0 +1,596 @@ +//////////////////////////////////////////////////////////////////////// +/** + * @file icaruscode/Timing/PMTBeamSignalsExtractor_module.cc + * @brief `icarus::timing::PMTBeamSignalsExtractor` producer module. + * @author Matteo Vicenzi (mvicenzi@bnl.gov) + * @date Sun Feb 11 11:37:14 2024 + **/ +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "art_root_io/TFileService.h" + +#include "sbnobj/Common/Trigger/ExtraTriggerInfo.h" +#include "icaruscode/Decode/ChannelMapping/IICARUSChannelMap.h" +#include "icaruscode/IcarusObj/PMTWaveformTimeCorrection.h" +#include "icaruscode/IcarusObj/PMTBeamSignal.h" +#include "lardataalg/DetectorInfo/DetectorTimingTypes.h" // electronics_time +#include "lardataalg/Utilities/quantities/spacetime.h" +#include "lardataobj/RawData/OpDetWaveform.h" +#include "lardataobj/RawData/TriggerData.h" + +#include "TTree.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace icarus::timing +{ + class PMTBeamSignalsExtractor; +} + +/** + * @brief Extracts RWM and EW times from special waveforms. + * + * This module extracts the RWM and EW from the waveforms recorded in the + * spare PMT channels and associates these times with all the channels + * in the same PMT readout crate. + * + * Input parameters + * ------------------------ + * + * This modules requires the following products: + * + * * `TriggerLabel` (input tag): tag for the `sbn::ExtraTriggerInfo` + * data product. This is required to check if it's a beam event or not. + * * `RWMlabel` (input tag): tag for the `std::vector` + * product containing the RWM special waveforms. + * * `EWlabel` (input tag): tag for the `std::vector` + * product containign the EW special waveforms. + * * `TriggerCorrectionLabel` (input tag): tag for the + * std::vector product that + * stores the channel-by-channel waveform timing corrections + * * `BoardSetup` (fhicl parameter set): description of the current + * V1730 setup from `CAEN_V1730_setup_icarus.fcl` mapping the special signals. + * It is meant to be the same configuration as used by the PMT decoding + * (see `BoardSetup` configuration parameter in `icarus::DaqDecoderICARUSPMT`). + + * * `ADCThreshold` (int): detection threshold to avoid cross-talk + * noise if one signal is missing from its waveform. + * * `DebugTrees` (bool): flag to produce plain ROOT trees for debugging. + * * `SaveWaveforms` (bool): flag to save full waveforms in the debug trees. + * + * Signal timing extraction + * ------------------------- + * + * The algorithm is quite unsophisticated and follows what is done for + * the digitized trigger signal in `PMTWaveformTimeCorrectionExtractor`. + * The reference signal is expected to be a sharp square wave in negative + * polarity. The time is based on the front side of that wave: + * + * * the absolute minimum of the waveform is found + * * an interval starting 20 ticks before that minimum is considered + * * the baseline level is defined as the value at the start of that interval + * * if the baseline-minimum difference is below a threshold, it is assumed to be noise + * and no time is returned + * * if not, the start time is set to the exact tick with an amplitude exceeding 20% + * of the minimum of the signal from the baseline + * + * + * Output products + * ------------------------- + * + * This modules produces two `std::vector` + * with 360 elements each, representing the relevent RWM or EW time for + * the corresponding PMT channel. When the discrimination algorithm decides + * the "peak" on a waveform to be noise, the corresponding entries are placed + * into the vector with an invalid value (`isValid()` returns `false`, + * but the identification data fields are correctly set). + * + * If the event is offbeam, these vectors are produced empty. + * + * * + * Debugging tree + * --------------- + * + * There are two debugging trees, enabled by `DebugTrees`, one for the Early Warning (EW) signals, + * named after the instance name of `EWlabel`, and one for the Resistive Wall Monitor (RWM) signals, + * named after the instance name of `RWMlabel`. Each entry in the tree represents a single signal + * in an event (thus, typically for ICARUS there are 8 entries per event). + * The content of an entry includes: + * * `run`, `event`: identifier of the event. + * * `timestamp`: timestamp of the event (in UTC), truncated at the second. + * * `n_channels`: the number of special signal channels in the event. + * * `channel`: the special channel nummber assigned to this signal. + * * `wfstart`: (uncorrected) waveform time start. + * * `sample`: the sample number in the waveform at which the signal was found. + * * `utime_abs`: uncorrected signal time [µs] + * * `time_abs`: corrected signal time [µs] + * * `time`: corrected signal time (relative to the trigger time) [µs] + * * `nsize` (only if `SaveWaveforms` is set): number of samples in waveform + * * `wf` (only if `SaveWaveforms` is set): full waveform. + */ + +class icarus::timing::PMTBeamSignalsExtractor : public art::EDProducer +{ +public: + explicit PMTBeamSignalsExtractor(fhicl::ParameterSet const &pset); + + // process waveforms + void extractBeamSignalTime(art::Event &e, art::InputTag const &label); + template + static T Median(std::vector data); + template + static std::size_t getMaxBin(std::vector const &vv, std::size_t startElement, std::size_t endElement); + template + static std::size_t getMinBin(std::vector const &vv, std::size_t startElement, std::size_t endElement); + template + static std::size_t getStartSample(std::vector const &vv, T thres); + + // unpack the V1730 special channels settings in a useful way + static std::map extractBoardBySpecialChannel(std::vector const &setup); + + // associate times to PMT channels + void associateBeamSignalsToChannels(art::InputTag const &label); + + // trigger-hardware timing correction + double getTriggerCorrection(int channel) const; + + // quick mapping conversions + std::string getDigitizerLabel(int channel) const; + std::string getCrate(int channel) const; + + // Plugins should not be copied or assigned. + PMTBeamSignalsExtractor(PMTBeamSignalsExtractor const &) = delete; + PMTBeamSignalsExtractor(PMTBeamSignalsExtractor &&) = delete; + PMTBeamSignalsExtractor &operator=(PMTBeamSignalsExtractor const &) = delete; + PMTBeamSignalsExtractor &operator=(PMTBeamSignalsExtractor &&) = delete; + + void beginJob() override; + void beginRun(art::Run &run) override; + void produce(art::Event &e) override; + +private: + /// Channel mappping + icarusDB::IICARUSChannelMap const &fChannelMap; + /// Save plain ROOT TTrees for debugging + bool const fDebugTrees; + /// Save raw waveforms in debug TTrees + bool const fSaveWaveforms; + /// Trigger instance label + art::InputTag const fTriggerLabel; + /// RWM waveform instance label + art::InputTag const fRWMlabel; + /// EW waveform instance label + art::InputTag const fEWlabel; + /// Trigger-hardware correction instance + art::InputTag const fTriggerCorrectionLabel; + /// Threshold for pulse selection + short int const fADCThreshold; + /// Special channel to board association in a map + std::map const fBoardBySpecialChannel; + + /// PMT sample duration [µs] + static constexpr double fPMTsamplingTick = 0.002; + /// Number of PMT channels + static constexpr std::size_t fNPMTChannels = 360; + + std::vector fCorrections; + std::map fBoardEffFragmentID; + double ftriggerTime; + + // output TTrees + std::map fOutTree; + // event info + int m_run; + int m_event; + int m_timestamp; + // special signal info + int m_n_channels; + unsigned int m_channel; + double m_wfstart; + std::size_t m_sample; + double m_time; + double m_time_abs; + double m_utime_abs; + std::size_t m_nsize; + std::vector m_wf; + + // prepare pointers for data products + using BeamSignalCollection = std::vector; + using BeamSignalCollectionPtr = std::unique_ptr; + + std::map> fBeamSignals; + std::map fSignalCollection; +}; + +// ----------------------------------------------------------------------------- + +icarus::timing::PMTBeamSignalsExtractor::PMTBeamSignalsExtractor(fhicl::ParameterSet const &pset) + : EDProducer{pset}, + fChannelMap(*(art::ServiceHandle{})), + fDebugTrees(pset.get("DebugTrees")), + fSaveWaveforms(pset.get("SaveWaveforms")), + fTriggerLabel(pset.get("TriggerLabel")), + fRWMlabel(pset.get("RWMlabel")), + fEWlabel(pset.get("EWlabel")), + fTriggerCorrectionLabel(pset.get("TriggerCorrectionLabel")), + fADCThreshold(pset.get("ADCThreshold")), + fBoardBySpecialChannel(extractBoardBySpecialChannel(pset.get>("BoardSetup"))) +{ + // Call appropriate consumes<>() functions here. + consumes>(fEWlabel); + consumes>(fRWMlabel); + + // Call appropriate produces<>() functions here. + produces>("RWM"); + produces>("EW"); +} + +// ----------------------------------------------------------------------------- + +void icarus::timing::PMTBeamSignalsExtractor::beginJob() +{ + // prepare outupt TTrees if requested + if (!fDebugTrees) + return; + + art::ServiceHandle tfs; + std::array const labels{fRWMlabel.instance(), fEWlabel.instance()}; + for (auto l : labels) + { + std::string name = l + "tree"; + std::string desc = l + " info"; + TTree *tree = tfs->make(name.c_str(), desc.c_str()); + tree->Branch("run", &m_run); + tree->Branch("event", &m_event); + tree->Branch("timestamp", &m_timestamp); + tree->Branch("n_channels", &m_n_channels); + tree->Branch("channel", &m_channel); + tree->Branch("wfstart", &m_wfstart); + tree->Branch("sample", &m_sample); + tree->Branch("utime_abs", &m_utime_abs); + tree->Branch("time_abs", &m_time_abs); + tree->Branch("time", &m_time); + + // add std::vector with full waveforms + // this can quickly make the TTrees quite heavy + if (fSaveWaveforms) + { + tree->Branch("nsize", &m_nsize); + tree->Branch("wf", &m_wf); + } + fOutTree[l] = tree; + } +} + +// ----------------------------------------------------------------------------- + +void icarus::timing::PMTBeamSignalsExtractor::beginRun(art::Run &run) +{ + // pre-save the association between digitizer_label and effective fragment ID + // needs to be done at the begin of each run in case mapping changed + fBoardEffFragmentID.clear(); + + for (unsigned int fragid = 0; fragid < fChannelMap.nPMTfragmentIDs(); fragid++) + { + auto const &pmtinfo = fChannelMap.getPMTchannelInfo(fragid)[0]; // pick first pmt on board + fBoardEffFragmentID[pmtinfo.digitizerLabel] = fragid; + } +} + +// ----------------------------------------------------------------------------- + +void icarus::timing::PMTBeamSignalsExtractor::produce(art::Event &e) +{ + + // initialize the data products + fBeamSignals.clear(); + fSignalCollection[fRWMlabel.instance()] = std::make_unique(); + fSignalCollection[fEWlabel.instance()] = std::make_unique(); + + // extract meta event information + m_run = e.id().run(); + m_event = e.id().event(); + m_timestamp = e.time().timeHigh(); // precision to the second + ftriggerTime = e.getProduct>(fTriggerLabel).at(0).TriggerTime(); // us + + // this module should run on beam-only events + // check the current beam gate + auto const &triggerInfo = e.getProduct(fTriggerLabel); + sbn::triggerSource const gateType = triggerInfo.sourceType; + + switch (gateType) + { + case sbn::triggerSource::BNB: + break; + case sbn::triggerSource::NuMI: + break; + default: + mf::LogTrace("PMTBeamSignalsExtractor") << "Skipping offbeam gate '" << name(gateType) << "'"; + e.put(std::move(fSignalCollection[fRWMlabel.instance()]), "RWM"); + e.put(std::move(fSignalCollection[fEWlabel.instance()]), "EW"); + return; + } + + // reside the collections to match the expected 360 PMT channels + // this is to avoid dynamic resizing later on + for (auto &collection : fSignalCollection) + collection.second->resize(fNPMTChannels); + + // get the trigger-hardware corrections that are applied on all signal waveforms + // if EW or RWM is to be compared to the PMT signals, it must be applied on them as well + // it also takes care of board-to-board offsets (see SBN-doc-34631, slide 5) + // Note: this is a vector of 360 elements, one correction for each signal channel + fCorrections = e.getProduct>(fTriggerCorrectionLabel); + std::size_t ntrig = fCorrections.size(); + + if (ntrig < 1) + mf::LogError("PMTBeamSignalsExtractor") << "Not found PMTWaveformTimeCorrections with label '" + << fTriggerCorrectionLabel.instance() << "'"; + else if (ntrig < fNPMTChannels) + mf::LogError("PMTBeamSignalsExtractor") << "Missing " << fNPMTChannels - ntrig << " PMTWaveformTimeCorrections with label '" + << fTriggerCorrectionLabel.instance() << "'"; + + // now the main course: getting EW and RWM waveforms + // information is stored by PMT crate name + extractBeamSignalTime(e, fRWMlabel); + extractBeamSignalTime(e, fEWlabel); + + // associating the proper RWM and EW time to each PMT channel + // collections are vectors of 360 elements (one value for each channel) + associateBeamSignalsToChannels(fRWMlabel); + associateBeamSignalsToChannels(fEWlabel); + + // place data products in the stream + // fix the cable swap for part of Run 2 right here!! + // see SBN-doc-34631 for details + if (gateType == sbn::triggerSource::BNB && m_run > 9704 && m_run < 11443) + std::swap(fSignalCollection[fRWMlabel.instance()], fSignalCollection[fEWlabel.instance()]); + + e.put(std::move(fSignalCollection[fRWMlabel.instance()]), "RWM"); + e.put(std::move(fSignalCollection[fEWlabel.instance()]), "EW"); +} + +// ----------------------------------------------------------------------------- + +void icarus::timing::PMTBeamSignalsExtractor::extractBeamSignalTime(art::Event &e, art::InputTag const &label) +{ + + std::string const &l = label.instance(); + auto const &waveforms = e.getProduct>(label); + m_n_channels = waveforms.size(); + + if (m_n_channels < 1) + mf::LogError("PMTBeamSignalsExtractor") << "Not found raw::OpDetWaveform with label '" << label.encode() << "'"; + else if (m_n_channels < 8) + mf::LogError("PMTBeamSignalsExtractor") << "Missing " << 8 - m_n_channels << " raw::OpDetWaveform with label '" + << label.encode() << "'"; + + m_wf.clear(); + + // get the start sample of the signals, one instance per PMT crate + // use threshold to skip spikes from electric crosstalk see SBN-doc-34928, slides 4-5. + // if no signal is found, set both times to icarus::timing::NoTime + // so that `isValid()` in PMTBeamSignal will complain if someone checks it + for (auto const &wave : waveforms) + { + + detinfo::timescales::electronics_time tstart = util::quantities::points::microsecond{wave.TimeStamp()}; + + // if nothing is found, first sample is returned (0) + m_sample = getStartSample(wave.Waveform(), fADCThreshold); + + m_channel = wave.ChannelNumber(); + m_wfstart = tstart.value(); + m_utime_abs = (m_sample != icarus::timing::NoSample) ? tstart.value() + fPMTsamplingTick * m_sample : icarus::timing::NoTime; + m_time_abs = (m_sample != icarus::timing::NoSample) ? m_utime_abs + getTriggerCorrection(m_channel) : icarus::timing::NoTime; + m_time = (m_sample != icarus::timing::NoSample) ? m_time_abs - ftriggerTime : icarus::timing::NoTime; + + std::string crate = getCrate(m_channel); + icarus::timing::PMTBeamSignal beamTime{m_channel, getDigitizerLabel(m_channel), crate, m_sample, m_time_abs, m_time}; + fBeamSignals[l].emplace(crate, std::move(beamTime)); + + if (fDebugTrees) + { + if (fSaveWaveforms) + { + m_wf = wave.Waveform(); + m_nsize = m_wf.size(); + } + fOutTree[l]->Fill(); + } + } +} + +// --------------------------------------------------------------------------- + +template +T icarus::timing::PMTBeamSignalsExtractor::Median(std::vector data) +{ + + std::nth_element(data.begin(), data.begin() + data.size() / 2, data.end()); + return data[data.size() / 2]; +} + +// ----------------------------------------------------------------------------- + +template +std::size_t icarus::timing::PMTBeamSignalsExtractor::getMinBin( + std::vector const &vv, std::size_t startElement, std::size_t endElement) +{ + + auto minel = + std::min_element(vv.begin() + startElement, vv.begin() + endElement); + std::size_t minsample = std::distance(vv.begin() + startElement, minel); + + return minsample; +} + +// ----------------------------------------------------------------------------- + +template +std::size_t icarus::timing::PMTBeamSignalsExtractor::getMaxBin( + std::vector const &vv, std::size_t startElement, std::size_t endElement) +{ + + auto maxel = + std::max_element(vv.begin() + startElement, vv.begin() + endElement); + + std::size_t maxsample = std::distance(vv.begin() + startElement, maxel); + + return maxsample; +} + +// ----------------------------------------------------------------------------- + +template +std::size_t icarus::timing::PMTBeamSignalsExtractor::getStartSample(std::vector const &vv, T thres) +{ + + // We are thinking in inverted polarity + std::size_t minbin = getMinBin(vv, 0, vv.size()); + + // Search only a cropped region of the waveform backward from the min + std::size_t maxbin = (minbin - 20) ? (minbin - 20) : 0; + + // Now we crawl betweem maxbin and minbin and we stop when: + // maxbin value - bin value > (maxbin value - minbin value )*0.2 + std::size_t startbin = maxbin; + auto delta = vv[maxbin] - vv[minbin]; + + if (delta < thres) // just noise + return icarus::timing::NoSample; // return no sample + + for (std::size_t bin = maxbin; bin < minbin; bin++) + { + auto val = vv[maxbin] - vv[bin]; + if (val >= 0.2 * delta) + { // 20% + startbin = bin - 1; + break; + } + } + + if (startbin < maxbin) + { + startbin = maxbin; + } + + return startbin; +} + +// ----------------------------------------------------------------------------- + +std::map icarus::timing::PMTBeamSignalsExtractor::extractBoardBySpecialChannel( + std::vector const &setup) +{ + + // map from special PMT channel to corresponding board + std::map boardBySpecialChannel; + + // unpack the V1730 special channels settings in a useful way + for (fhicl::ParameterSet const &s : setup) + { + auto const &innerSet = s.get>("SpecialChannels"); + boardBySpecialChannel[innerSet[0].get("Channel")] = s.get("Name"); + } + + return boardBySpecialChannel; +} + +// ----------------------------------------------------------------------------- + +std::string icarus::timing::PMTBeamSignalsExtractor::getDigitizerLabel(int channel) const +{ + + // get the board name, convert to digitizer_label + std::string board = fBoardBySpecialChannel.at(channel); // eg. icaruspmtewtop02 + + std::string head = "icaruspmt"; + char letter = board.back() - '1' + 'A'; // converts 01,02,03 to A,B,C + + board.erase(board.find(head), head.size()); // eg. ewtop02 + std::transform(board.begin(), board.end(), board.begin(), ::toupper); // eg. EWTOP02 + board.insert(2, 1, '-'); // insert dash at position 2, e.g: EW-TOP02 + board.insert(6, 1, '-'); // insert dash at position 6, e.g: EW-TOP-02 + + return board.substr(0, board.size() - 2) + letter; // e.g: EW-TOP-B +} + +// ----------------------------------------------------------------------------- + +std::string icarus::timing::PMTBeamSignalsExtractor::getCrate(int channel) const +{ + + std::string digitizer_label = getDigitizerLabel(channel); + return digitizer_label.substr(0, digitizer_label.size() - 2); +} + +// ----------------------------------------------------------------------------- + +double icarus::timing::PMTBeamSignalsExtractor::getTriggerCorrection(int channel) const +{ + + std::string digitizer_label = getDigitizerLabel(channel); + + // trigger-hardware corrections are shared by all channels on the same board + // we can pick the first channel on the desired board + int fragID = fBoardEffFragmentID.at(digitizer_label); + auto pmtinfo = fChannelMap.getPMTchannelInfo(fragID)[0]; // pick first ch on board + int pmtch = pmtinfo.channelID; + + // trigger-hardware correction are in order + // index of vector is pmtch + return fCorrections.at(pmtch).startTime; +} + +// ----------------------------------------------------------------------------- + +void icarus::timing::PMTBeamSignalsExtractor::associateBeamSignalsToChannels(art::InputTag const &label) +{ + + std::string const &l = label.instance(); + + // loop through the signals which are one per PMT crate + // for each crate, find the corresponding digitizers + for (auto const &signal : fBeamSignals[l]) + { + + // build the PMT digitizer labels that live in this crate + // then convert it into fragment id + std::array const letters{"-A", "-B", "-C"}; + for (auto letter : letters) + { + + std::string digitizer_label = signal.first + letter; + int fragID = fBoardEffFragmentID[digitizer_label]; + + // use the fragment id to access the PMT channels + // mapping works via fragment id (it's very annoying..) + // loop through the PMT channels and set their RWM/EW time + for (auto pmtinfo : fChannelMap.getPMTchannelInfo(fragID)) + { + + std::size_t channel = pmtinfo.channelID; + // collections are resized to fNPMTChannels + // always room in the collection (channel -> vector index) + fSignalCollection[l]->at(channel) = signal.second; + + } // for each channel + } // for each board + } // for each crate +} + +DEFINE_ART_MODULE(icarus::timing::PMTBeamSignalsExtractor) diff --git a/icaruscode/Timing/timing_beam.fcl b/icaruscode/Timing/timing_beam.fcl new file mode 100644 index 000000000..a68836ff5 --- /dev/null +++ b/icaruscode/Timing/timing_beam.fcl @@ -0,0 +1,18 @@ +#include "CAEN_V1730_setup_icarus.fcl" + +BEGIN_PROLOG + +icarus_beam_signal_extractor: +{ + module_type: "PMTBeamSignalsExtractor" + TriggerLabel: "daqTrigger" + EWlabel: "daqPMT:EW" + RWMlabel: "daqPMT:RWM" + TriggerCorrectionLabel: "daqPMT:globtrg" + ADCThreshold: 200 + BoardSetup: @local::icarus_V1730_setup + DebugTrees: false + SaveWaveforms: false +} + +END_PROLOG