From ab4671453554b27a19557c097f8d7f7d70668fb2 Mon Sep 17 00:00:00 2001 From: jzennamo Date: Thu, 1 Oct 2026 13:20:31 -0500 Subject: [PATCH 1/2] fixes for beam accounting --- sbncode/BeamSpillInfoRetriever/CMakeLists.txt | 2 +- .../ICARUSBNBRetriever_module.cc | 80 ++-- .../BeamSpillInfoRetriever/MWRMatching.cpp | 43 +++ sbncode/BeamSpillInfoRetriever/MWRMatching.h | 50 +++ sbncode/BeamSpillInfoRetriever/POTTools.cpp | 97 +++-- sbncode/BeamSpillInfoRetriever/POTTools.h | 25 +- .../SBNDBNBRetriever_module.cc | 68 ++-- .../SBNDBNBZEROBIASRetriever_module.cc | 65 ++-- sbncode/BeamSpillInfoRetriever/getFOM.cpp | 362 +++++++----------- sbncode/BeamSpillInfoRetriever/getFOM.h | 17 +- .../job/icarusbnbspillinfo.fcl | 2 + .../job/sbndbnbdefaults.fcl | 2 + sbncode/CAFMaker/FillExposure.cxx | 8 +- 13 files changed, 422 insertions(+), 399 deletions(-) create mode 100644 sbncode/BeamSpillInfoRetriever/MWRMatching.cpp create mode 100644 sbncode/BeamSpillInfoRetriever/MWRMatching.h diff --git a/sbncode/BeamSpillInfoRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt index 5332b931d..38be9044b 100644 --- a/sbncode/BeamSpillInfoRetriever/CMakeLists.txt +++ b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt @@ -34,7 +34,7 @@ art_make_library( larcorealg::CoreUtils LIBRARY_NAME sbn_POTTools - SOURCE POTTools.cpp + SOURCE POTTools.cpp MWRMatching.cpp ) art_make_library( diff --git a/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc index cee5fa7a5..dee476e90 100644 --- a/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc @@ -86,6 +86,16 @@ class sbn::ICARUSBNBRetriever : public art::EDProducer { Comment{ "Window of time in seconds to use for mwr ifbeam queries." } }; + fhicl::Atom MWRMaxTimeDiff { + Name{ "MWRMaxTimeDiff" }, + Comment{ "largest time difference between a multiwire reading and its spill [s]; <= 0 disables" }, + 0.0333 // default: half the 15 Hz Booster period + }; + fhicl::Atom ReuseLastBPMOffsets { + Name{ "ReuseLastBPMOffsets" }, + Comment{ "if a BPM offset cannot be read, use the last valid value seen in this job" }, + true // default + }; fhicl::Atom TriggerDatabaseFile { Name{ "TriggerDatabaseFile" }, Comment{ "path of local database of all recorded events and their trigger, in SQLite format" } @@ -135,7 +145,10 @@ class sbn::ICARUSBNBRetriever : public art::EDProducer { sqlite3 *db; int rc; - static constexpr double MWRtoroidDelay = -0.035; ///< the same time point is measured _t_ by MWR and _t + MWRtoroidDelay`_ by the toroid [ms] + static constexpr double MWRtoroidDelay = -0.035; ///< the same time point is measured _t_ by MWR and _t + MWRtoroidDelay`_ by the toroid [s] + double fMWRMaxTimeDiff; ///< largest accepted multiwire-spill time difference [s] + bool fReuseLastBPMOffsets; ///< whether to fill missing BPM offsets from fOffsetCache + mutable sbn::pot::BPMOffsetCache_t fOffsetCache; ///< last valid BPM offsets seen /// Returns the information of the trigger in the current event. sbn::pot::TriggerInfo_t extractTriggerInfo(art::Event const& e) const; @@ -229,6 +242,8 @@ sbn::ICARUSBNBRetriever::ICARUSBNBRetriever(Parameters const& params) vp873( ifbeam_handle->getBeamFolder(params().VP873Bundle(), params().URL(), params().TimeWindow())), offsets( ifbeam_handle->getBeamFolder(params().OffsetBundle(), params().URL(), params().TimeWindow())), bfp_mwr( ifbeam_handle->getBeamFolder(params().MultiWireBundle(), params().URL(), params().MWR_TimeWindow())), + fMWRMaxTimeDiff(params().MWRMaxTimeDiff()), + fReuseLastBPMOffsets(params().ReuseLastBPMOffsets()), fTriggerDatabaseFile(params().TriggerDatabaseFile()) { @@ -366,8 +381,6 @@ int sbn::ICARUSBNBRetriever::matchMultiWireData( // DAQ trigger times int spill_count = 0; int spills_removed = 0; - std::vector matched_MWR; - matched_MWR.resize(3); ///reject time_stamps which have a trigger_type == 1 from data-base //To-Do @@ -410,47 +423,14 @@ int sbn::ICARUSBNBRetriever::matchMultiWireData( //Great we found a matched spill! Let's count it spill_count++; - //Loop through the multiwire devices: - - for(int dev = 0; dev < int(MWR_times.size()); dev++){ - - //Loop through the multiwire times: - double Tdiff = 1000000000.; - matched_MWR[dev] = 0; - - for(int mwrt = 0; mwrt < int(MWR_times[dev].size()); mwrt++){ - - //found a candidate match! - if(fabs((MWR_times[dev][mwrt] - times_temps[i])) >= Tdiff){continue;} - - bool best_match = true; - - //Check for a better match... - for (size_t j = 0; j < times_temps.size(); j++) { - if( j == i) continue; - if(times_temps[j] > (triggerInfo.t_current_event+fTimePad)){continue;} - if(times_temps[j] <= (triggerInfo.t_previous_event+fTimePad)){continue;} - - //is there a better match later in the spill sequence - if(fabs((MWR_times[dev][mwrt] - times_temps[j])) < - fabs((MWR_times[dev][mwrt] - times_temps[i]))){ - //we can have patience... - best_match = false; - break; - } - }//end better match check - - //Verified best match! - if(best_match == true){ - matched_MWR[dev] = mwrt; - Tdiff = fabs((MWR_times[dev][mwrt] - times_temps[i])); - } - - }//end loop over MWR times - - }//end loop over MWR devices - - sbn::BNBSpillInfo spillInfo = sbn::pot::makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873); + // Associate one reading of each multiwire device to this spill (-1: none) + std::vector const matched_MWR = sbn::pot::matchMWRToSpill( + MWR_times, times_temps, i, + triggerInfo.t_previous_event+fTimePad, triggerInfo.t_current_event+fTimePad, + fMWRMaxTimeDiff); + + sbn::BNBSpillInfo spillInfo = sbn::pot::makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873, + fReuseLastBPMOffsets? &fOffsetCache: nullptr); std::tuple allFOM = sbn::getBNBqualityFOM(spillInfo); spillInfo.FOM = std::get<0>(allFOM); spillInfo.PreFitFOM = std::get<1>(allFOM); @@ -484,6 +464,18 @@ void sbn::ICARUSBNBRetriever::endSubRun(art::SubRun& sr) mf::LogDebug("ICARUSBNBRetriever")<< "Total number of DAQ Spills : " << TotalBeamSpills << std::endl; mf::LogDebug("ICARUSBNBRetriever")<< "Total number of Selected Spills : " << fOutbeamInfos.size() << std::endl; + // Spills recorded before the first valid BPM offset reading of the job could + // not be patched from the cache at the time: do it now and redo their FOM. + if (fReuseLastBPMOffsets) { + for (sbn::BNBSpillInfo& info: fOutbeamInfos) { + if (!sbn::pot::fillMissingBPMOffsets(info, fOffsetCache)) continue; + std::tuple const allFOM = sbn::getBNBqualityFOM(info); + info.FOM = std::get<0>(allFOM); + info.PreFitFOM = std::get<1>(allFOM); + info.NoMultiWireFOM = std::get<2>(allFOM); + } + } + auto p = std::make_unique< std::vector< sbn::BNBSpillInfo > >(); std::swap(*p, fOutbeamInfos); diff --git a/sbncode/BeamSpillInfoRetriever/MWRMatching.cpp b/sbncode/BeamSpillInfoRetriever/MWRMatching.cpp new file mode 100644 index 000000000..1aa3cff5a --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/MWRMatching.cpp @@ -0,0 +1,43 @@ +/** + * @file sbncode/BeamSpillInfoRetriever/MWRMatching.cpp + * @brief Association of multiwire chamber readings to BNB spills. + */ +#include "sbncode/BeamSpillInfoRetriever/MWRMatching.h" + +#include +#include + +std::vector sbn::pot::matchMWRToSpill(std::vector> const& MWR_times, + std::vector const& spill_times, + std::size_t i, + double windowLow, double windowHigh, + double maxTimeDiff) +{ + std::vector matched(MWR_times.size(), -1); + if (i >= spill_times.size()) return matched; + double const t_spill = spill_times[i]; + + for (std::size_t dev = 0; dev < MWR_times.size(); ++dev) { + double Tdiff = std::numeric_limits::max(); + for (std::size_t mwrt = 0; mwrt < MWR_times[dev].size(); ++mwrt) { + double const t_mwr = MWR_times[dev][mwrt]; + double const d = std::abs(t_mwr - t_spill); + if (d >= Tdiff) continue; + + // is another spill in the window a better match for this reading? + bool best_match = true; + for (std::size_t j = 0; j < spill_times.size(); ++j) { + if (j == i) continue; + if (spill_times[j] > windowHigh) continue; + if (spill_times[j] <= windowLow) continue; + if (std::abs(t_mwr - spill_times[j]) < d) { best_match = false; break; } + } + if (best_match) { + matched[dev] = static_cast(mwrt); + Tdiff = d; + } + } + if (matched[dev] >= 0 && maxTimeDiff > 0. && Tdiff > maxTimeDiff) matched[dev] = -1; + } + return matched; +} diff --git a/sbncode/BeamSpillInfoRetriever/MWRMatching.h b/sbncode/BeamSpillInfoRetriever/MWRMatching.h new file mode 100644 index 000000000..8648b6e3f --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/MWRMatching.h @@ -0,0 +1,50 @@ +#ifndef SBNCODE_BEAMSPILLINFORETRIEVER_MWRMATCHING_H +#define SBNCODE_BEAMSPILLINFORETRIEVER_MWRMATCHING_H + +/** + * @file sbncode/BeamSpillInfoRetriever/MWRMatching.h + * @brief Association of multiwire chamber readings to BNB spills. + * + * Factored out of the ICARUS/SBND BNB retriever modules, which carried three + * copies of the same loop. Kept free of art dependencies so it can be tested + * standalone. + */ + +#include +#include + +namespace sbn::pot { + + /** + * @brief Finds, for each multiwire device, the reading belonging to spill `i`. + * @param MWR_times reading times per device [s] (already corrected for the + * MWR-to-toroid delay) + * @param spill_times toroid spill times [s] + * @param i index in `spill_times` of the spill to match + * @param windowLow spills with time <= windowLow are not considered as + * competitors + * @param windowHigh spills with time > windowHigh are not considered as + * competitors + * @param maxTimeDiff largest accepted |t_MWR - t_spill| [s]; non-positive + * disables the requirement + * @return one index per device into `MWR_times[dev]`, or `-1` if no reading + * can be associated with this spill + * + * A reading is a candidate for spill `i` only if no other spill in the + * window is closer to it in time; among candidates the closest one wins. + * + * Differences from the original inline code: + * - no reading -> `-1` (it used to silently fall back to index 0, attaching + * an arbitrary, possibly far away, profile to the spill); + * - the best candidate is rejected if it is more than `maxTimeDiff` away + * (there was no limit: profiles tens of seconds away were being used). + */ + std::vector matchMWRToSpill(std::vector> const& MWR_times, + std::vector const& spill_times, + std::size_t i, + double windowLow, double windowHigh, + double maxTimeDiff); + +} // namespace sbn::pot + +#endif diff --git a/sbncode/BeamSpillInfoRetriever/POTTools.cpp b/sbncode/BeamSpillInfoRetriever/POTTools.cpp index a83d5f304..625e1630e 100644 --- a/sbncode/BeamSpillInfoRetriever/POTTools.cpp +++ b/sbncode/BeamSpillInfoRetriever/POTTools.cpp @@ -89,7 +89,7 @@ namespace sbn::pot{ } sbn::BNBSpillInfo makeBNBSpillInfo - (art::EventID const& eventID, double time, MWRdata_t const& MWRdata, std::vector const& matched_MWR, std::unique_ptr const& bfp, const std::unique_ptr & offsets,const std::unique_ptr & vp873) + (art::EventID const& eventID, double time, MWRdata_t const& MWRdata, std::vector const& matched_MWR, std::unique_ptr const& bfp, const std::unique_ptr & offsets,const std::unique_ptr & vp873, BPMOffsetCache_t* offsetCache) { auto const& [ MWR_times, unpacked_MWR ] = MWRdata; // alias @@ -164,13 +164,27 @@ namespace sbn::pot{ try{bfp->GetNamedData(time, "E:M876HM",&M876HM);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} try{bfp->GetNamedData(time, "E:M876VM",&M876VM);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} - try{offsets->GetNamedData(time, "E_VP873S",&VP873Offset);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} - try{offsets->GetNamedData(time, "E_HP875S",&HP875Offset);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} - try{offsets->GetNamedData(time, "E_VP875S",&VP875Offset);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} - try{offsets->GetNamedData(time, "E_HPTG1S",&HPTG1Offset);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} - try{offsets->GetNamedData(time, "E_HPTG2S",&HPTG2Offset);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} - try{offsets->GetNamedData(time, "E_VPTG1S",&VPTG1Offset);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} - try{offsets->GetNamedData(time, "E_VPTG2S",&VPTG2Offset);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + // BPM offsets are settings that change rarely; the query (epsilon of a few + // minutes) does fail at times, which used to leave -999 and cost the spill + // its FOM. With a cache, fall back to the last valid value. + auto getOffset = [&](std::string const& name, double& value) { + try{offsets->GetNamedData(time, name.c_str(), &value);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + if (!offsetCache) return; + if (value != -999) { + (*offsetCache)[name] = value; + } + else if (auto const it = offsetCache->find(name); it != offsetCache->end()) { + value = it->second; + mf::LogDebug("BNBRetriever") << "At time : " << time << " " << name << " missing, using last valid value " << value; + } + }; + getOffset("E_VP873S", VP873Offset); + getOffset("E_HP875S", HP875Offset); + getOffset("E_VP875S", VP875Offset); + getOffset("E_HPTG1S", HPTG1Offset); + getOffset("E_HPTG2S", HPTG2Offset); + getOffset("E_VPTG1S", VPTG1Offset); + getOffset("E_VPTG2S", VPTG2Offset); //crunch the times unsigned long int time_closest_int = (int) TOR860_time; @@ -214,37 +228,25 @@ namespace sbn::pot{ beamInfo.FOM = FOM; - for(auto const& MWRdata: unpacked_MWR){ - std::ignore = MWRdata; - assert(!MWRdata.empty()); - } - - if(unpacked_MWR[0].empty()){ - beamInfo.M875BB.clear(); - beamInfo.M875BB_spill_time_diff = -999;//units in seconds - } - else{ - beamInfo.M875BB = unpacked_MWR[0][matched_MWR[0]]; - beamInfo.M875BB_spill_time_diff = (MWR_times[0][matched_MWR[0]] - time); - } - - if(unpacked_MWR[1].empty()){ - beamInfo.M876BB.clear(); - beamInfo.M876BB_spill_time_diff = -999;//units in seconds - } - else{ - beamInfo.M876BB = unpacked_MWR[1][matched_MWR[1]]; - beamInfo.M876BB_spill_time_diff = (MWR_times[1][matched_MWR[1]] - time); - } - - if(unpacked_MWR[2].empty()){ - beamInfo.MMBTBB.clear(); - beamInfo.MMBTBB_spill_time_diff = -999;//units in seconds - } - else{ - beamInfo.MMBTBB = unpacked_MWR[2][matched_MWR[2]]; - beamInfo.MMBTBB_spill_time_diff = (MWR_times[2][matched_MWR[2]] - time); - } + // A device has no reading for this spill if it reported nothing at all or + // if no reading could be associated with this spill (matched index -1). + // (The previous `assert(!MWRdata.empty())` on every device contradicted the + // handling of empty devices below and would abort debug builds.) + auto fillMWR = [&](std::size_t dev, auto& profile, auto& timeDiff) { + int const idx = (dev < matched_MWR.size())? matched_MWR[dev]: -1; + if (unpacked_MWR[dev].empty() || idx < 0 || std::size_t(idx) >= unpacked_MWR[dev].size()) { + profile.clear(); + timeDiff = -999; //units in seconds + } + else { + profile = unpacked_MWR[dev][idx]; + timeDiff = MWR_times[dev][idx] - time; + } + }; + fillMWR(0, beamInfo.M875BB, beamInfo.M875BB_spill_time_diff); + fillMWR(1, beamInfo.M876BB, beamInfo.M876BB_spill_time_diff); + fillMWR(2, beamInfo.MMBTBB, beamInfo.MMBTBB_spill_time_diff); + // We do not write these to the art::Events because // we can filter events but want to keep all the POT // information, so we'll write it to the SubRun @@ -258,6 +260,23 @@ namespace sbn::pot{ return beamInfo; } + bool fillMissingBPMOffsets(sbn::BNBSpillInfo& info, BPMOffsetCache_t const& cache) + { + bool filled = false; + auto fill = [&](char const* name, auto& value) { + if (value != -999) return; + if (auto const it = cache.find(name); it != cache.end()) { value = it->second; filled = true; } + }; + fill("E_VP873S", info.VP873Offset); + fill("E_HP875S", info.HP875Offset); + fill("E_VP875S", info.VP875Offset); + fill("E_HPTG1S", info.HPTG1Offset); + fill("E_HPTG2S", info.HPTG2Offset); + fill("E_VPTG1S", info.VPTG1Offset); + fill("E_VPTG2S", info.VPTG2Offset); + return filled; + } + bool BrokenClock(double time, std::unique_ptr const& bfp) { double TOR860 = -999; // units e12 protons diff --git a/sbncode/BeamSpillInfoRetriever/POTTools.h b/sbncode/BeamSpillInfoRetriever/POTTools.h index 1e54647a8..b2c3b4d27 100644 --- a/sbncode/BeamSpillInfoRetriever/POTTools.h +++ b/sbncode/BeamSpillInfoRetriever/POTTools.h @@ -12,7 +12,10 @@ #include "artdaq-core/Data/ContainerFragment.hh" #include "sbncode/BeamSpillInfoRetriever/MWRData.h" +#include "sbncode/BeamSpillInfoRetriever/MWRMatching.h" #include "larcorealg/CoreUtils/counter.h" +#include +#include #include namespace sbn::pot{ @@ -38,6 +41,9 @@ namespace sbn::pot{ std::vector< std::vector< std::vector< int > > > unpacked_MWR; } MWRdata_t; + /// Last valid value of each BPM offset (device name -> offset [mm]). + using BPMOffsetCache_t = std::map; + /** * @brief Extracts information from PTB for a single HLT for use in SBND POT accounting. * @@ -85,8 +91,25 @@ namespace sbn::pot{ MWRdata_t extractSpillTimes(TriggerInfo_t const& triggerInfo, std::unique_ptr const& bfp, std::unique_ptr const& bfp_mwr, double fTimePad, double MWRtoroidDelay, sbn::MWRData mwrdata ); /** * @brief Compile spill information into BNBSpillInfo object + * @param matched_MWR index of the multiwire reading of each device for this + * spill (see `matchMWRToSpill()`); `-1` means none + * @param offsetCache if not null, BPM offsets that cannot be read from the + * database are taken from the last valid reading (they are + * slowly changing settings), and valid readings update it + * + * Without `offsetCache` a failed offset query leaves the offset at -999, and + * the spill gets no FOM even when all the BPM readings are present. + */ + sbn::BNBSpillInfo makeBNBSpillInfo(art::EventID const& eventID, double time, MWRdata_t const& MWRdata, std::vector const& matched_MWR, std::unique_ptr const& bfp, std::unique_ptr const& offsets, std::unique_ptr const& vp873, BPMOffsetCache_t* offsetCache = nullptr); + + /** + * @brief Fills BPM offsets still missing (-999) in `info` from `cache`. + * @return whether any offset was filled (the FOM then needs recomputing) + * + * Meant for the end of the subrun: spills recorded before the first valid + * offset reading of the job could not use the cache when they were made. */ - sbn::BNBSpillInfo makeBNBSpillInfo(art::EventID const& eventID, double time, MWRdata_t const& MWRdata, std::vector const& matched_MWR, std::unique_ptr const& bfp, std::unique_ptr const& offsets, std::unique_ptr const& vp873); + bool fillMissingBPMOffsets(sbn::BNBSpillInfo& info, BPMOffsetCache_t const& cache); } #endif diff --git a/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc index b76da2059..8830f9cb0 100644 --- a/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc @@ -42,7 +42,10 @@ class sbn::SBNDBNBRetriever : public art::EDProducer { sbn::MWRData mwrdata; art::ServiceHandle ifbeam_handle; - static constexpr double MWRtoroidDelay = -0.035; ///< the same time point is measured _t_ by MWR and _t + MWRtoroidDelay`_ by the toroid [ms] + static constexpr double MWRtoroidDelay = -0.035; ///< the same time point is measured _t_ by MWR and _t + MWRtoroidDelay`_ by the toroid [s] + double fMWRMaxTimeDiff; ///< largest accepted multiwire-spill time difference [s] + bool fReuseLastBPMOffsets; ///< whether to fill missing BPM offsets from fOffsetCache + mutable sbn::pot::BPMOffsetCache_t fOffsetCache; ///< last valid BPM offsets seen sbn::pot::TriggerInfo_t extractTriggerInfo(art::Event const& e) const; int matchMultiWireData( @@ -60,6 +63,8 @@ sbn::SBNDBNBRetriever::SBNDBNBRetriever(fhicl::ParameterSet const & params) fTimePad = params.get("TimePadding"); fBESOffset = params.get("BESOffset"); fDeviceUsedForTiming = params.get("DeviceUsedForTiming"); + fMWRMaxTimeDiff = params.get("MWRMaxTimeDiff", 0.0333); + fReuseLastBPMOffsets = params.get("ReuseLastBPMOffsets", true); double const timeWindow = std::stod(params.get("TimeWindow")); bfp = ifbeam_handle->getBeamFolder(params.get("Bundle"), params.get("URL"), timeWindow); bfp->set_epsilon(0.02); @@ -163,8 +168,6 @@ int sbn::SBNDBNBRetriever::matchMultiWireData( // DAQ trigger times int spill_count = 0; int spills_removed = 0; - std::vector matched_MWR; - matched_MWR.resize(3); // Iterating through each of the beamline times for (size_t i = 0; i < times_temps.size(); i++) { @@ -181,47 +184,14 @@ int sbn::SBNDBNBRetriever::matchMultiWireData( //Great we found a matched spill! Let's count it spill_count++; - //Loop through the multiwire devices: + // Associate one reading of each multiwire device to this spill (-1: none) + std::vector const matched_MWR = sbn::pot::matchMWRToSpill( + MWR_times, times_temps, i, + triggerInfo.t_previous_event+fTimePad, triggerInfo.t_current_event+fTimePad, + fMWRMaxTimeDiff); - for(int dev = 0; dev < int(MWR_times.size()); dev++){ - - //Loop through the multiwire times: - double Tdiff = 1000000000.; - matched_MWR[dev] = 0; - - for(int mwrt = 0; mwrt < int(MWR_times[dev].size()); mwrt++){ - - //found a candidate match! - if(fabs((MWR_times[dev][mwrt] - times_temps[i])) >= Tdiff){continue;} - - bool best_match = true; - - //Check for a better match... - for (size_t j = 0; j < times_temps.size(); j++) { - if( j == i) continue; - if(times_temps[j] > (triggerInfo.t_current_event+fTimePad)){continue;} - if(times_temps[j] <= (triggerInfo.t_previous_event+fTimePad)){continue;} - - //is there a better match later in the spill sequence - if(fabs((MWR_times[dev][mwrt] - times_temps[j])) < - fabs((MWR_times[dev][mwrt] - times_temps[i]))){ - //we can have patience... - best_match = false; - break; - } - }//end better match check - - //Verified best match! - if(best_match == true){ - matched_MWR[dev] = mwrt; - Tdiff = fabs((MWR_times[dev][mwrt] - times_temps[i])); - } - - }//end loop over MWR times - - }//end loop over MWR devices - - sbn::BNBSpillInfo spillInfo = makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873); + sbn::BNBSpillInfo spillInfo = makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873, + fReuseLastBPMOffsets? &fOffsetCache: nullptr); std::tuple allFOM = sbn::getBNBqualityFOM(spillInfo); spillInfo.FOM = std::get<0>(allFOM); spillInfo.PreFitFOM = std::get<1>(allFOM); @@ -248,6 +218,18 @@ void sbn::SBNDBNBRetriever::endSubRun(art::SubRun& sr) mf::LogDebug("SBNDBNBRetriever")<< "Total number of DAQ Spills : " << TotalBeamSpills << std::endl; mf::LogDebug("SBNDBNBRetriever")<< "Total number of Selected Spills : " << fOutbeamInfos.size() << std::endl; + // Spills recorded before the first valid BPM offset reading of the job could + // not be patched from the cache at the time: do it now and redo their FOM. + if (fReuseLastBPMOffsets) { + for (sbn::BNBSpillInfo& info: fOutbeamInfos) { + if (!sbn::pot::fillMissingBPMOffsets(info, fOffsetCache)) continue; + std::tuple const allFOM = sbn::getBNBqualityFOM(info); + info.FOM = std::get<0>(allFOM); + info.PreFitFOM = std::get<1>(allFOM); + info.NoMultiWireFOM = std::get<2>(allFOM); + } + } + auto p = std::make_unique< std::vector< sbn::BNBSpillInfo > >(); std::swap(*p, fOutbeamInfos); diff --git a/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc index f301eff14..d0e88bc25 100644 --- a/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc @@ -40,7 +40,10 @@ class sbn::SBNDBNBZEROBIASRetriever : public art::EDProducer { std::unique_ptr bfp_mwr; sbn::MWRData mwrdata; art::ServiceHandle ifbeam_handle; - static constexpr double MWRtoroidDelay = -0.035; ///< the same time point is measured _t_ by MWR and _t + MWRtoroidDelay`_ by the toroid [ms] + static constexpr double MWRtoroidDelay = -0.035; ///< the same time point is measured _t_ by MWR and _t + MWRtoroidDelay`_ by the toroid [s] + double fMWRMaxTimeDiff; ///< largest accepted multiwire-spill time difference [s] + bool fReuseLastBPMOffsets; ///< whether to fill missing BPM offsets from fOffsetCache + mutable sbn::pot::BPMOffsetCache_t fOffsetCache; ///< last valid BPM offsets seen std::vector< sbn::BNBSpillInfo > fOutbeamInfos; std::vector< sbn::BNBSpillInfo > fOutbeamInfosTotal; @@ -59,6 +62,8 @@ sbn::SBNDBNBZEROBIASRetriever::SBNDBNBZEROBIASRetriever(fhicl::ParameterSet cons fTimePad = params.get("TimePadding"); fDeviceUsedForTiming = params.get("DeviceUsedForTiming"); fBESOffset = params.get("BESOffset"); + fMWRMaxTimeDiff = params.get("MWRMaxTimeDiff", 0.0333); + fReuseLastBPMOffsets = params.get("ReuseLastBPMOffsets", true); const double timeWindow = std::stod(params.get("TimeWindow")); bfp = ifbeam_handle->getBeamFolder(params.get("Bundle"), params.get("URL"), timeWindow); bfp->set_epsilon(0.02); @@ -165,14 +170,15 @@ void sbn::SBNDBNBZEROBIASRetriever::matchMultiWireData( // We'll keep track of how many of these spills match to our // DAQ trigger times int spills_removed = 0; - std::vector matched_MWR; - matched_MWR.resize(3); - // Iterating through each of the beamline times - + // Pick the latest spill in the window (closest before the current trigger). + // Fixed: the broken-clock check tested times_temps[i] (the previous best, + // initially entry 0) instead of the candidate times_temps[k]; and when no + // spill qualified, entry 0 (out of range for an empty list) was used anyway. double best_diff = 10000000000.0; double diff; size_t i = 0; + bool found = false; for (size_t k = 0; k < times_temps.size(); k++){ diff = (triggerInfo.t_current_event + fTimePad) - times_temps[k];// diff is greater than zero! @@ -184,50 +190,29 @@ void sbn::SBNDBNBZEROBIASRetriever::matchMultiWireData( spills_removed++; continue;} - if(sbn::pot::BrokenClock(times_temps[i], bfp)){ + if(sbn::pot::BrokenClock(times_temps[k], bfp)){ continue; } best_diff = diff; i = k; + found = true; } } - for(int dev = 0; dev < int(MWR_times.size()); dev++){ - //Loop through the multiwire times: - double Tdiff = 1000000000.; - matched_MWR[dev] = 0; - - for(int mwrt = 0; mwrt < int(MWR_times[dev].size()); mwrt++){ - //found a candidate match! - if(fabs((MWR_times[dev][mwrt] - times_temps[i])) >= Tdiff){continue;} - - bool best_match = true; - - for (size_t j = 0; j < times_temps.size(); j++) { - //Check for a better match... - if( j == i) continue; - if(times_temps[j] > (triggerInfo.t_current_event+fTimePad)){continue;} - if(times_temps[j] <= (triggerInfo.t_previous_event+fTimePad)){continue;} - - //is there a better match later in the spill sequence - if(fabs((MWR_times[dev][mwrt] - times_temps[j])) < - fabs((MWR_times[dev][mwrt] - times_temps[i]))){ - //we can have patience... - best_match = false; - break; - } - }//end better match check - - //Verified best match! - if(best_match == true){ - matched_MWR[dev] = mwrt; - Tdiff = fabs((MWR_times[dev][mwrt] - times_temps[i])); - } - }//end loop over MWR times - }//end loop over MWR devices + if (!found) { + mf::LogDebug("SBNDBNBZEROBIASRetriever") << "matchMultiWireData:: no spill found in the window, skipping event " << eventID; + return; + } + + // Associate one reading of each multiwire device to this spill (-1: none) + std::vector const matched_MWR = sbn::pot::matchMWRToSpill( + MWR_times, times_temps, i, + triggerInfo.t_previous_event+fTimePad, triggerInfo.t_current_event+fTimePad, + fMWRMaxTimeDiff); - sbn::BNBSpillInfo spillInfo = makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873); + sbn::BNBSpillInfo spillInfo = makeBNBSpillInfo(eventID, times_temps[i], MWRdata, matched_MWR, bfp, offsets, vp873, + fReuseLastBPMOffsets? &fOffsetCache: nullptr); std::tuple allFOM = sbn::getBNBqualityFOM(spillInfo); spillInfo.FOM = std::get<0>(allFOM); spillInfo.PreFitFOM = std::get<1>(allFOM); diff --git a/sbncode/BeamSpillInfoRetriever/getFOM.cpp b/sbncode/BeamSpillInfoRetriever/getFOM.cpp index e8447cfae..c7bebfd90 100644 --- a/sbncode/BeamSpillInfoRetriever/getFOM.cpp +++ b/sbncode/BeamSpillInfoRetriever/getFOM.cpp @@ -5,6 +5,9 @@ */ #include "sbncode/BeamSpillInfoRetriever/getFOM.h" #include +#include +#include +#include #include "TH1D.h" #include "TFitResult.h" #include @@ -15,32 +18,38 @@ using namespace std; namespace sbn { - bool onePlot = true; + namespace { - std::tuple getBNBqualityFOM(BNBSpillInfo & spill ) - { - double fom=0; - double prefitfom=0; - double noMWfom=0; + /// Value used by the retrievers for a device that could not be read. + constexpr double MissingValue = -999.; + + /// A device reading is usable if it is finite and not the "missing" marker. + bool isValid(double v) { return std::isfinite(v) && v != MissingValue; } + + /// Number of wires per plane in the multiwire chambers. + constexpr std::size_t NWires = 48; + + /// Straight-line extrapolation to the target centre through two BPMs. + /// Positions in mm, z in m: the slope is in mm/m, i.e. mrad, which is what + /// calcFOM() expects for the angle (no atan: the slope already *is* the angle). + std::pair extrapolate + (double delta0, double z0, double delta1, double z1, double ztarget) + { + double const ang = (delta1 - delta0) / (z1 - z0); + double const pos = delta0 + ang * (ztarget - z0); + return { ang, pos }; + } + + } // local namespace - double hp875_offset= spill.HP875Offset; - double vp875_offset= spill.VP875Offset; - double vp873_offset= spill.VP873Offset; - double hptg1_offset= spill.HPTG1Offset; - //double vptg1_offset= spill.VPTG1Offset; - double hptg2_offset= spill.HPTG2Offset; - double vptg2_offset= spill.VPTG2Offset; - - //Decides which position monitor to check first - int useHTG = 1; - int useVTG = 1; - //Z Position of the monitors in m + std::tuple getBNBqualityFOM(BNBSpillInfo const& spill) + { + //Z Position of the monitors in m double const vp873_zpos= 191.153656; - double const hp875_zpos= 202.116104; + double const hp875_zpos= 202.116104; double const vp875_zpos= 202.3193205; double const hptg1_zpos= 204.833267; - //double const vptg1_zpos= 204.629608; double const hptg2_zpos= 205.240662; double const vptg2_zpos= 205.036835; double const target_center_zpos= 206.870895; @@ -48,201 +57,102 @@ namespace sbn { double const p875y[]={0.279128, 0.337048, 0}; double const p876x[]={0.166172, 0.30999, -0.00630299}; double const p876y[]={0.13425, 0.580862, 0}; - - - std::vector tor860; - std::vector tor875; - std::vector hp875; - std::vector vp875; - std::vector vp873; - std::vector hptg1; - std::vector vptg1; - std::vector hptg2; - std::vector vptg2; - std::vector m875hs; // Multiwire station after Mag 875, Fit to Horizontal Sigma - std::vector m875hm; // Multiwire station after Mag 875, Fit to Horizontal Mean - std::vector m875vs; // Multiwire station after Mag 875, Fit to Vertical Sigma - std::vector m875vm; // Multiwire station after Mag 875, Fit to Vertical Mean + // ---- intensity: TOR860, falling back to TOR875. + // A missing toroid is stored as -999e12; a spill with no (or negative) + // intensity has no meaningful FOM (the optics model needs ppp > 0). + double tor = -1.; + if (isValid(spill.TOR860) && spill.TOR860 > 0.) tor = spill.TOR860; + else if (isValid(spill.TOR875) && spill.TOR875 > 0.) tor = spill.TOR875; + else return {-1, -1, -1}; - std::vector m876hs; // Multiwire station after Mag 876, Fit to Horizontal Sigma - std::vector m876hm; // Multiwire station after Mag 876, Fit to Horizontal Mean - std::vector m876vs; // Multiwire station after Mag 876, Fit to Vertical Sigma - std::vector m876vm; // Multiwire station after Mag 876, Fit to Vertical Mean - - std::vector mw875(spill.M875BB.begin(), spill.M875BB.end()); - std::vector mw876(spill.M876BB.begin(), spill.M876BB.end()); - std::vector mwtgt(spill.MMBTBB.begin(), spill.MMBTBB.end()); - - - tor860.push_back(spill.TOR860); - tor875.push_back(spill.TOR875); - hp875.push_back(spill.HP875); - vp875.push_back(spill.VP875); - vp873.push_back(spill.VP873); - hptg1.push_back(spill.HPTG1); - vptg1.push_back(spill.VPTG1); - hptg2.push_back(spill.HPTG2); - vptg2.push_back(spill.VPTG2); - - m875hs.push_back(spill.M875HS); - m875hm.push_back(spill.M875HM); - m875vs.push_back(spill.M875VS); - m875vm.push_back(spill.M875VM); - - m876hs.push_back(spill.M876HS); - m876hm.push_back(spill.M876HM); - m876vs.push_back(spill.M876VS); - m876vm.push_back(spill.M876VM); - - double tor; - if (!tor860.empty()) - tor=tor860[0]; - else if (!tor875.empty()) - tor=tor875[0]; - else - return {-1,-1,-1}; - - /** - * @brief when creating ntuples for pot counting script the variables are filled with -999 - * this could create a difference when passing events with FOM>1 since - * events with missing BPM data would get FOM=2, while events with BPM set to -999 will get FOM=0 - * bad or missing MWR data gets FOM 4 in either case - */ - if (hptg2.empty()) hptg2.push_back(-999); - if (hptg1.empty()) hptg1.push_back(-999); - if (hp875.empty()) hp875.push_back(-999); - if (vptg2.empty()) vptg2.push_back(-999); - if (vptg1.empty()) vptg1.push_back(-999); - if (vp875.empty()) vp875.push_back(-999); - if (vp873.empty()) vp873.push_back(-999); - double horang; - - auto interpolate_hp875 = [delta_hp875=(hp875[0]-hp875_offset), hp875_zpos, target_center_zpos] - (double delta_value, double zpos) - { - double const ang = (delta_value-delta_hp875)/(zpos-hp875_zpos); - double const pos = delta_hp875+ang*(target_center_zpos-hp875_zpos); - return std::pair(ang, pos); - }; - - // return 2 when missing essential beam horizontal position data: - if (hp875.empty() || (hptg1.empty() && hptg2.empty())) return {2,2,2}; - bool const doUseHTG1 = (useHTG == 1) || hptg2.empty(); - auto const [ Tanhorang, horpos ] = doUseHTG1? - interpolate_hp875(hptg1[0] - hptg1_offset, hptg1_zpos): - interpolate_hp875(hptg2[0] - hptg2_offset, hptg2_zpos); + // ---- horizontal: HP875 and HPTG1 (HPTG2 as fallback), offsets subtracted. + // Missing devices used to be fed into the extrapolation as -999 (the + // previous `.empty()` checks could never trigger), placing the beam ~1 m + // off target and losing the spill. + bool const okHP875 = isValid(spill.HP875) && isValid(spill.HP875Offset); + bool const okHPTG1 = isValid(spill.HPTG1) && isValid(spill.HPTG1Offset); + bool const okHPTG2 = isValid(spill.HPTG2) && isValid(spill.HPTG2Offset); + if (!okHP875 || (!okHPTG1 && !okHPTG2)) return {2, 2, 2}; + double const delta_hp875 = spill.HP875 - spill.HP875Offset; + auto const [ horang, horpos ] = okHPTG1 + ? extrapolate(delta_hp875, hp875_zpos, spill.HPTG1 - spill.HPTG1Offset, hptg1_zpos, target_center_zpos) + : extrapolate(delta_hp875, hp875_zpos, spill.HPTG2 - spill.HPTG2Offset, hptg2_zpos, target_center_zpos); + // ---- vertical: VP875 and VP873 (VPTG2 as fallback), offsets subtracted. + bool const okVP875 = isValid(spill.VP875) && isValid(spill.VP875Offset); + bool const okVP873 = isValid(spill.VP873) && isValid(spill.VP873Offset); + bool const okVPTG2 = isValid(spill.VPTG2) && isValid(spill.VPTG2Offset); + if (!okVP875 || (!okVP873 && !okVPTG2)) return {3, 3, 3}; + double const delta_vp875 = spill.VP875 - spill.VP875Offset; + auto const [ verang, verpos ] = okVP873 + ? extrapolate(delta_vp875, vp875_zpos, spill.VP873 - spill.VP873Offset, vp873_zpos, target_center_zpos) + : extrapolate(delta_vp875, vp875_zpos, spill.VPTG2 - spill.VPTG2Offset, vptg2_zpos, target_center_zpos); - double verang; - auto interpolate_vp875 = [delta_vp875=(vp875[0]-vp875_offset), vp875_zpos, target_center_zpos] - (double delta_value, double zpos) - { - double const ang = (delta_value-delta_vp875)/(zpos-vp875_zpos); - double const pos = delta_vp875+ang*(target_center_zpos-vp875_zpos); - return std::pair(ang, pos); - }; - - // return 3 when missing essential beam horizontal position data: - if (vp875.empty() || (vptg1.empty() && vptg2.empty())) return {3,3,3}; - bool const doUseVTG1 = (useVTG == 1) || vptg2.empty(); - auto const [ Tanverang, verpos ] = doUseVTG1? - interpolate_vp875(vp873[0] - vp873_offset, vp873_zpos): - interpolate_vp875(vptg2[0] - vptg2_offset, vptg2_zpos); + const double smallSigmaX =0.5, largeSigmaX = 10, smallSigmaY = 0.3, largeSigmaY =10, maxChi2X = 20, maxChi2Y = 20; + auto inWindow = [&](double sx, double sy) + { return sx>smallSigmaX && sxsmallSigmaY && sy0) { - processBNBprofile(&mwtgt[FirstXMWtgt], xx, sx,chi2x); - processBNBprofile(&mwtgt[FirstYMWtgt], yy, sy, chi2y); - if (sx>smallSigmaX && sxsmallSigmaY && sy 0`, + // allowed reading past the end of a short vector) and both fits succeed. + // Preference: target multiwire (no transformation), then M876, then M875. + struct MWDevice_t { + std::vector const* data; + double const* px; double const* py; + }; + MWDevice_t const mwDevices[] = { + { &spill.MMBTBB, nullptr, nullptr }, + { &spill.M876BB, p876x, p876y }, + { &spill.M875BB, p875x, p875y }, + }; + double tgtsx = MissingValue, tgtsy = MissingValue; + bool goodFit = false; + for (auto const& dev: mwDevices) { + if (dev.data->size() < 2*NWires) continue; + std::vector const mw(dev.data->begin(), dev.data->begin() + 2*NWires); + double xx, yy, sx, sy, chi2x, chi2y; + bool const fitOK = processBNBprofile(&mw[0], xx, sx, chi2x) + & processBNBprofile(&mw[NWires], yy, sy, chi2y); + if (!fitOK) continue; + if (dev.px) { + sx = dev.px[0] + dev.px[1]*sx + dev.px[2]*sx*sx; + sy = dev.py[0] + dev.py[1]*sy + dev.py[2]*sy*sy; } - } - if (!good_tgt && mw876.size()>0) { - processBNBprofile(&mw876[FirstXMWtgt], xx,sx,chi2x); - processBNBprofile(&mw876[FirstYMWtgt], yy,sy,chi2y); - double tgtsx876=p876x[0]+p876x[1]*sx+p876x[2]*sx*sx; - double tgtsy876=p876y[0]+p876y[1]*sy+p876y[2]*sy*sy; - if (tgtsx876>smallSigmaX && tgtsx876smallSigmaY && tgtsy8760){ - processBNBprofile(&mw875[FirstXMWtgt], xx,sx,chi2x); - processBNBprofile(&mw875[FirstYMWtgt], yy,sy,chi2y); - double tgtsx875=p875x[0]+p875x[1]*sx+p875x[2]*sx*sx; - double tgtsy875=p875y[0]+p875y[1]*sy+p875y[2]*sy*sy; - if (tgtsx875>smallSigmaX && tgtsx875smallSigmaY && tgtsy875smallSigmaX && tgtsx876smallSigmaY && tgtsy876smallSigmaX && tgtsx875smallSigmaY && tgtsy875SetBinContent(i+1,-mwdata[i]-minx); + for (unsigned int i=0;i threshold && first_x==-1) first_x=i; if (-mwdata[i]-minx > threshold) last_x=i+1; - hProf->SetBinError(i+1,error); - } - if (hProf->GetSumOfWeights()>0) { - TFitResultPtr const fit = hProf->Fit("gaus","QNS","",-12+first_x*0.5,-12+last_x*0.5); - x = fit->Parameter(1); - sx = fit->Parameter(2); - chi2= fit->Chi2() / fit->Ndf(); - delete hProf; - } else { - x=99999; - sx=99999; - chi2=99999; + hProf.SetBinError(i+1,error); } - entry += 1; + if (hProf.GetSumOfWeights() <= 0 || first_x < 0) return false; + + TFitResultPtr const fit = hProf.Fit("gaus","QNS","",-12+first_x*0.5,-12+last_x*0.5); + if (!fit.Get() || fit->Status() != 0 || fit->Ndf() <= 0) return false; + x = fit->Parameter(1); + sx = fit->Parameter(2); + chi2= fit->Chi2() / fit->Ndf(); + return true; } - - + + double calcFOM(double horpos, double horang, double verpos, double verang, double ppp, double tgtsx, double tgtsy) { ppp /= 1e12; //converts to 10^12 POT @@ -480,7 +390,7 @@ namespace sbn { } x = x + dx; } - sum = sum*dx*dy/(2.0*3.14159*sx*sy*sqrt(1.0-rho2)); + sum = sum*dx*dy/(2.0*M_PI*sx*sy*sqrt(1.0-rho2)); // add a guard for double precision diff --git a/sbncode/BeamSpillInfoRetriever/getFOM.h b/sbncode/BeamSpillInfoRetriever/getFOM.h index 78fd0de22..a3cba4ac9 100644 --- a/sbncode/BeamSpillInfoRetriever/getFOM.h +++ b/sbncode/BeamSpillInfoRetriever/getFOM.h @@ -9,6 +9,8 @@ #include "sbnobj/Common/POTAccounting/BNBSpillInfo.h" +#include + namespace sbn @@ -18,8 +20,18 @@ namespace sbn * * The figure of merit is described in [SBN DocDB 41901](https://sbn-docdb.fnal.gov/cgi-bin/sso/ShowDocument?docid=41901). * Inputs the BNBSpillInfo and returns the BNB Quality Metric called FOM, derived from MicroBooNE's FOM + * + * @return { FOM with fitted multiwire width, FOM with database ("pre-fit") + * width, FOM with nominal width } + * + * Each value is in [0, 1], or one of these codes (same for all three): + * * `-1`: no valid intensity (TOR860 and TOR875 missing or <= 0, i.e. no beam) + * * `2`: horizontal position/angle cannot be formed (missing BPM or BPM offset) + * * `3`: vertical position/angle cannot be formed (missing BPM or BPM offset) + * The first two values are additionally `-999` when no usable width is found. + * A device value of `-999` is treated as missing. */ - std::tuple getBNBqualityFOM(BNBSpillInfo& spill); + std::tuple getBNBqualityFOM(BNBSpillInfo const& spill); /** * @brief Inside the getFOM script, takes the positions and angles of the beam and calculates the BNB FOM @@ -41,8 +53,9 @@ namespace sbn /** * @brief Inputs the MWR Data and determines the centroid, sigma, and chi2 value of a gaussian fit of the beam + * @return whether the fit succeeded (valid status, positive degrees of freedom) */ - void processBNBprofile(const double* mwdata, double &x, double& sx, double& chi2); + bool processBNBprofile(const double* mwdata, double &x, double& sx, double& chi2); } #endif diff --git a/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl b/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl index 03c63ab3e..30d39e36c 100644 --- a/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl +++ b/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl @@ -16,5 +16,7 @@ icarusbnbspillinfo: { raw_data_label: "daq" DeviceUsedForTiming: "E:TOR860" TriggerDatabaseFile: "triggerDatabase/icarus_triggers.db" + MWRMaxTimeDiff: 0.0333 #unit seconds, multiwire readings further than this from the spill are not used (<= 0 disables) + ReuseLastBPMOffsets: true #if a BPM offset query fails, use the last valid offset seen in this job } END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl b/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl index 47f37d9ba..524586264 100644 --- a/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl +++ b/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl @@ -14,6 +14,8 @@ sbndbnbspillinfo: { TimeWindow: "700" #seconds MWR_TimeWindow: "700" #seconds DeviceUsedForTiming: "E:TOR860" + MWRMaxTimeDiff: 0.0333 #unit seconds, multiwire readings further than this from the spill are not used (<= 0 disables) + ReuseLastBPMOffsets: true #if a BPM offset query fails, use the last valid offset seen in this job } END_PROLOG diff --git a/sbncode/CAFMaker/FillExposure.cxx b/sbncode/CAFMaker/FillExposure.cxx index 91746f26b..12decbfa8 100644 --- a/sbncode/CAFMaker/FillExposure.cxx +++ b/sbncode/CAFMaker/FillExposure.cxx @@ -14,15 +14,17 @@ namespace caf std::cout << "makeSRBNBInfo: Manually calculated width:" << FOM << std::endl; std::cout << "makeSRBNBInfo: Pre-Fit width from database:" << PreFitFOM << std::endl; std::cout << "makeSRBNBInfo: Assuming width of 1.0:" << NoMultiWireFOM << std::endl; - if((FOM > 0.0) & (FOM <= 1.0)){ + // valid FOM values are in [0, 1]; getBNBqualityFOM() flags failures with + // -999, -1, 2, 3 (a FOM of exactly 0, beam fully off target, is valid) + if((FOM >= 0.0) && (FOM <= 1.0)){ finalFOM = FOM; std::cout << "makeSRBNBInfo: Chose manually calculated width" << std::endl; } - else if((PreFitFOM > 0.0) & (PreFitFOM <= 1.0)){ + else if((PreFitFOM >= 0.0) && (PreFitFOM <= 1.0)){ finalFOM = PreFitFOM; std::cout << "makeSRBNBInfo: Chose pre-fit width" << std::endl; } - else if((NoMultiWireFOM > 0.0) & (NoMultiWireFOM <= 1.0)) + else if((NoMultiWireFOM >= 0.0) && (NoMultiWireFOM <= 1.0)) { finalFOM = 100.+NoMultiWireFOM; std::cout << "makeSRBNBInfo: Chose assumed width of 1.0" << std::endl; From 32fdb95c0d28e69240093a96b10409a86f455883 Mon Sep 17 00:00:00 2001 From: jzennamo Date: Tue, 6 Oct 2026 10:39:53 -0500 Subject: [PATCH 2/2] adding the fixes --- sbncode/BeamSpillInfoRetriever/BNBFOMFill.cpp | 180 +++++++++++ sbncode/BeamSpillInfoRetriever/BNBFOMFill.h | 74 +++++ sbncode/BeamSpillInfoRetriever/CMakeLists.txt | 2 +- .../ICARUSBNBRetriever_module.cc | 29 ++ .../SBNDBNBRetriever_module.cc | 14 + .../SBNDBNBZEROBIASRetriever_module.cc | 14 + sbncode/BeamSpillInfoRetriever/getFOM.cpp | 306 ++++++++++++------ sbncode/BeamSpillInfoRetriever/getFOM.h | 55 ++++ .../job/icarusbnbspillinfo.fcl | 3 + .../job/sbndbnbdefaults.fcl | 3 + sbncode/CAFMaker/FillExposure.cxx | 7 +- 11 files changed, 582 insertions(+), 105 deletions(-) create mode 100644 sbncode/BeamSpillInfoRetriever/BNBFOMFill.cpp create mode 100644 sbncode/BeamSpillInfoRetriever/BNBFOMFill.h diff --git a/sbncode/BeamSpillInfoRetriever/BNBFOMFill.cpp b/sbncode/BeamSpillInfoRetriever/BNBFOMFill.cpp new file mode 100644 index 000000000..702ff2085 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/BNBFOMFill.cpp @@ -0,0 +1,180 @@ +/** + * @file sbncode/BeamSpillInfoRetriever/BNBFOMFill.cpp + * @brief BNB figure of merit for spills with missing inputs (see BNBFOMFill.h). + */ +#include "sbncode/BeamSpillInfoRetriever/BNBFOMFill.h" +#include "sbncode/BeamSpillInfoRetriever/getFOM.h" + +#include +#include +#include + +namespace sbn { + + namespace { + + constexpr double MissingValue = -999.; + bool isValid(double v) { return std::isfinite(v) && v != MissingValue; } + + double spillTime(BNBSpillInfo const& s) { return s.spill_time_s + 1e-9 * s.spill_time_ns; } + + bool measuredWidth(BNBBeamState const& s) { return s.widthSource != BNBBeamState::Nominal; } + + /// Sets the FOM fields of a spill that was filled or got a new width. + void storeFOM(BNBSpillInfo& spill, BNBBeamState const& state, double fom) { + if (measuredWidth(state)) { + spill.FOM = fom; + spill.PreFitFOM = MissingValue; + spill.NoMultiWireFOM = computeFOM(state, true); + } + else { + spill.FOM = spill.PreFitFOM = MissingValue; + spill.NoMultiWireFOM = fom; + } + } + + } // local namespace + + + std::vector improveBNBqualityFOMs + (std::vector& spills, BNBFOMFillConfig const& cfg) + { + using namespace fomstatus; + std::size_t const n = spills.size(); + + // spills in time order + std::vector order(n); + std::iota(order.begin(), order.end(), 0); + std::stable_sort(order.begin(), order.end(), + [&spills](std::size_t a, std::size_t b){ return spillTime(spills[a]) < spillTime(spills[b]); }); + + std::vector S(n); + std::vector fom(n, MissingValue), t(n); + for (std::size_t j = 0; j < n; ++j) { + BNBSpillInfo const& spill = spills[order[j]]; + S[j] = getBNBBeamState(spill); + if (S[j].hasFOM()) fom[j] = computeFOM(S[j]); + t[j] = spillTime(spill); + } + auto const valid = [&](std::size_t j){ return S[j].hasFOM(); }; + // the spills just before and after j, both within nbMaxGap + auto const adjacent = [&](std::size_t j){ + return j > 0 && j + 1 < n && t[j] - t[j-1] <= cfg.nbMaxGap && t[j+1] - t[j] <= cfg.nbMaxGap; + }; + std::vector changed(n, false); + + // ---- 1. width of an empty multiwire from the adjacent spills + if (cfg.neighbourWidth) { + auto const goodWidth = [&](std::size_t j) + { return valid(j) && measuredWidth(S[j]) && !(S[j].status & WidthFromNeighbors); }; + for (std::size_t j = 0; j < n; ++j) { + if (!valid(j) || !(S[j].status & MWEmpty) || !adjacent(j)) continue; + BNBBeamState const &p = S[j-1], &q = S[j+1]; + if (!goodWidth(j-1) || !goodWidth(j+1)) continue; + if (std::abs(p.sx - q.sx) > cfg.nbMaxDSig || std::abs(p.sy - q.sy) > cfg.nbMaxDSig) continue; + S[j].sx = 0.5 * (p.sx + q.sx); + S[j].sy = 0.5 * (p.sy + q.sy); + S[j].widthSource = p.widthSource; + S[j].status = (S[j].status & ~NoMWWidth) | WidthFromNeighbors; + fom[j] = computeFOM(S[j]); + changed[j] = true; + } + } + + // ---- 2. neighbour fill: a BPM reading is missing, the adjacent spills are stable + if (cfg.neighbourFill) { + std::vector const S0 = S; // neighbours as they were + std::vector const fom0 = fom; + std::vector filled(n, false); + for (std::size_t j = 0; j < n; ++j) { + BNBBeamState& s = S[j]; + if (valid(j) || (s.status & NoTOR) || s.tor <= cfg.minTor || !adjacent(j)) continue; + if (!(s.status & (NoHBPM | NoVBPM))) continue; + BNBBeamState const &p = S0[j-1], &q = S0[j+1]; + if (!p.hasFOM() || !q.hasFOM() || filled[j-1]) continue; + double const tm = 0.5 * (p.tor + q.tor); + if (std::abs(s.tor - tm) / tm > cfg.nbMaxDTor) continue; + if (std::max(std::abs(p.hpos - q.hpos), std::abs(p.vpos - q.vpos)) > cfg.nbMaxDPos) continue; + if (std::max(std::abs(p.hang - q.hang), std::abs(p.vang - q.vang)) > cfg.nbMaxDAng) continue; + if (measuredWidth(p) && measuredWidth(q) // a nominal width on either side: no requirement + && std::max(std::abs(p.sx - q.sx), std::abs(p.sy - q.sy)) > cfg.nbMaxDSig) continue; + if (std::abs(fom0[j-1] - fom0[j+1]) > cfg.nbMaxDFOM) continue; + if (cfg.nbMinFOM >= 0. && (fom0[j-1] <= cfg.nbMinFOM || fom0[j+1] <= cfg.nbMinFOM)) continue; + s.hpos = 0.5 * (p.hpos + q.hpos); s.hang = 0.5 * (p.hang + q.hang); + s.vpos = 0.5 * (p.vpos + q.vpos); s.vang = 0.5 * (p.vang + q.vang); + if (!measuredWidth(s)) { // no width of its own + if (measuredWidth(p) && measuredWidth(q)) { + s.sx = 0.5 * (p.sx + q.sx); s.sy = 0.5 * (p.sy + q.sy); + s.widthSource = p.widthSource; + s.status = (s.status & ~NoMWWidth) | WidthFromNeighbors; + } + } + // the "missing" bits stay set: they record why the spill was filled + s.status |= NeighborFilled; + BNBBeamState computed = s; + computed.status &= ~(NoHBPM | NoVBPM); + fom[j] = computeFOM(computed); + filled[j] = changed[j] = true; + } + } + + // ---- 3. 875-station drop-outs: own target BPMs + the measured spills on both sides + if (cfg.burstFill) { + auto const isGood = [&](std::size_t j){ + BNBSpillInfo const& sp = spills[order[j]]; + return valid(j) && !(S[j].status & NeighborFilled) + && isValid(sp.HP875) && isValid(sp.VP875) && isValid(sp.VP873) + && isValid(sp.HPTG1) && isValid(sp.VPTG2); + }; + std::vector good; + for (std::size_t j = 0; j < n; ++j) if (isGood(j)) good.push_back(j); + for (std::size_t j = 0; j < n && !good.empty(); ++j) { + BNBBeamState& s = S[j]; + BNBSpillInfo& sp = spills[order[j]]; + if (valid(j) || (s.status & (NoTOR | NeighborFilled)) || s.tor <= cfg.minTor) continue; + if (isValid(sp.HP875) || isValid(sp.VP875) || !isValid(sp.HPTG1) || !isValid(sp.VPTG2)) continue; + auto const it = std::lower_bound(good.begin(), good.end(), j); // first good after j + if (it == good.begin() || it == good.end()) continue; + std::size_t const a = *(it - 1), b = *it; + if (t[j] - t[a] > cfg.burstMaxDt || t[b] - t[j] > cfg.burstMaxDt) continue; + if (fom[a] <= cfg.burstMinFOM || fom[b] <= cfg.burstMinFOM) continue; + BNBBeamState const &A = S[a], &B = S[b]; + if (std::abs(A.hang - B.hang) >= cfg.burstMaxDAng || std::abs(A.vang - B.vang) >= cfg.burstMaxDAng) continue; + BNBSpillInfo const &spA = spills[order[a]], &spB = spills[order[b]]; + if (std::abs(sp.HPTG1 - 0.5 * (spA.HPTG1 + spB.HPTG1)) >= cfg.burstMaxDTgt) continue; + if (std::abs(sp.VPTG2 - 0.5 * (spA.VPTG2 + spB.VPTG2)) >= cfg.burstMaxDTgt) continue; + BNBBeamState e = s; + e.hpos = 0.5 * ((A.hpos - spA.HPTG1) + (B.hpos - spB.HPTG1)) + sp.HPTG1; + e.hang = 0.5 * (A.hang + B.hang); + e.vpos = 0.5 * ((A.vpos - spA.VPTG2) + (B.vpos - spB.VPTG2)) + sp.VPTG2; + e.vang = 0.5 * (A.vang + B.vang); + // own width unless the M876 database width says the chamber was empty + bool const m876empty = (isValid(sp.M876HS) && sp.M876HS > 4.0) || (isValid(sp.M876VS) && sp.M876VS > 4.0); + if (m876empty || !measuredWidth(e)) { + e.sx = e.sy = MissingValue; + e.widthSource = BNBBeamState::Nominal; + e.status |= NoMWWidth; + } + e.status &= ~(NoHBPM | NoVBPM); + double const f = computeFOM(e); + if (!(f > cfg.burstAccept)) continue; + unsigned int const missing = s.status & (NoHBPM | NoVBPM); // kept: why the spill was filled + s = e; + s.status |= BurstFill | missing; + fom[j] = f; + changed[j] = true; + } + } + + std::vector status(n); + for (std::size_t j = 0; j < n; ++j) { + status[order[j]] = S[j].status; + if (!changed[j]) continue; + BNBBeamState computed = S[j]; + computed.status &= ~(NoHBPM | NoVBPM); + storeFOM(spills[order[j]], computed, fom[j]); + } + return status; + } + +} // namespace sbn diff --git a/sbncode/BeamSpillInfoRetriever/BNBFOMFill.h b/sbncode/BeamSpillInfoRetriever/BNBFOMFill.h new file mode 100644 index 000000000..6649c015c --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/BNBFOMFill.h @@ -0,0 +1,74 @@ +#ifndef SBNCODE_BEAMSPILLRETRIEVER_BNBFOMFILL_H +#define SBNCODE_BEAMSPILLRETRIEVER_BNBFOMFILL_H + +/** + * @file sbncode/BeamSpillInfoRetriever/BNBFOMFill.h + * @brief BNB figure of merit for spills with missing inputs, from their neighbours. + * + * These use the other spills of the same subrun, so they run at the end of + * the subrun in the retriever modules. Validated with closure tests on + * ICARUS Run 2 and SBND Run 1 (hide what the fill replaces in measured spills, + * rebuild the FOM, count bad spills, true FOM <= 0.98, that get accepted). + */ + +#include "sbnobj/Common/POTAccounting/BNBSpillInfo.h" + +#include + +namespace sbn { + + struct BNBFOMFillConfig { + + double minTor = 1e11; ///< spills below this intensity [protons] are not filled + + // ---- empty multiwire / neighbour width + /// Spills whose database width is an empty chamber take the width of the + /// adjacent spills when those have a measured width agreeing within + /// `nbMaxDSig` (otherwise they keep the nominal width). + bool neighbourWidth = true; + + // ---- neighbour fill: BPM reading missing, adjacent spills stable + bool neighbourFill = true; + double nbMaxGap = 0.1; ///< each adjacent spill at most this far [s] (next 15 Hz spill) + double nbMaxDTor = 0.03; ///< spill intensity within this fraction of the neighbours' mean + double nbMaxDPos = 1.0; ///< neighbours' positions agree [mm] + double nbMaxDAng = 0.4; ///< neighbours' angles agree [mrad] + double nbMaxDSig = 0.1; ///< neighbours' widths agree [mm] + double nbMaxDFOM = 0.001; ///< neighbours' FOMs agree + double nbMinFOM = -1.; ///< both neighbours' FOM above this (< 0: no requirement; ICARUS: 0.998) + + // ---- 875-station drop-outs (HP875 and VP875 missing, target BPMs read) + bool burstFill = false; ///< ICARUS + double burstMaxDt = 60.; ///< measured spill on each side within this [s] + double burstMinFOM = 0.999; ///< both of them above this FOM + double burstMaxDAng = 0.2; ///< their angles agree [mrad] + double burstMaxDTgt = 0.5; ///< own HPTG1/VPTG2 within this of their mean [mm] + double burstAccept = 0.995; ///< filled only if the estimated FOM is above this + + }; // BNBFOMFillConfig + + /** + * @brief Improves the FOM of the spills of one subrun using their neighbours. + * @param spills all the spills of the subrun (any order); FOMs are updated + * @param config which fills to apply and their thresholds + * @return the status word (`sbn::fomstatus` bits) of each spill, same order + * + * Spills that are not filled keep the result of `getBNBqualityFOM()`. + * A filled spill gets its FOM in `FOM` (measured width) or in + * `NoMultiWireFOM` (nominal width), the others set to -999; its BPM + * readings are not changed, so a filled spill still shows missing readings. + * + * 1. neighbour width (bit `WidthFromNeighbors`) + * 2. neighbour fill (bit `NeighborFilled`): the spill takes the mean position + * and angle of the spills just before and after (and their width if it + * has none), with its own intensity + * 3. 875 drop-out fill (bit `BurstFill`): own HPTG1/VPTG2 plus the offset to + * the target and the angle of the nearest measured spill on each side; + * BPM offsets cancel in the differences + */ + std::vector improveBNBqualityFOMs + (std::vector& spills, BNBFOMFillConfig const& config); + +} // namespace sbn + +#endif // SBNCODE_BEAMSPILLRETRIEVER_BNBFOMFILL_H diff --git a/sbncode/BeamSpillInfoRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt index 38be9044b..99f531734 100644 --- a/sbncode/BeamSpillInfoRetriever/CMakeLists.txt +++ b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt @@ -46,7 +46,7 @@ art_make_library( LIBRARY_NAME sbn_getFOM - SOURCE getFOM.cpp + SOURCE getFOM.cpp BNBFOMFill.cpp ) install_headers() diff --git a/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc index dee476e90..356304b7f 100644 --- a/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/ICARUSBNBRetriever/ICARUSBNBRetriever_module.cc @@ -18,6 +18,7 @@ #include #include "sbncode/BeamSpillInfoRetriever/POTTools.h" #include "sbncode/BeamSpillInfoRetriever/getFOM.h" +#include "sbncode/BeamSpillInfoRetriever/BNBFOMFill.h" namespace sbn { class ICARUSBNBRetriever; @@ -96,6 +97,21 @@ class sbn::ICARUSBNBRetriever : public art::EDProducer { Comment{ "if a BPM offset cannot be read, use the last valid value seen in this job" }, true // default }; + fhicl::Atom ImproveFOM { + Name{ "ImproveFOM" }, + Comment{ "at the end of the subrun, recover the FOM of spills with missing inputs from their neighbours (see BNBFOMFill.h)" }, + true // default + }; + fhicl::Atom NeighbourMinFOM { + Name{ "NeighbourMinFOM" }, + Comment{ "neighbour fill only if both adjacent spills have FOM above this (< 0: no requirement)" }, + 0.998 // default (ICARUS Run 2 closure test) + }; + fhicl::Atom FillBPMDropouts { + Name{ "FillBPMDropouts" }, + Comment{ "fill HP875+VP875 drop-outs from the target BPMs and the measured spills on both sides" }, + true // default + }; fhicl::Atom TriggerDatabaseFile { Name{ "TriggerDatabaseFile" }, Comment{ "path of local database of all recorded events and their trigger, in SQLite format" } @@ -149,6 +165,8 @@ class sbn::ICARUSBNBRetriever : public art::EDProducer { double fMWRMaxTimeDiff; ///< largest accepted multiwire-spill time difference [s] bool fReuseLastBPMOffsets; ///< whether to fill missing BPM offsets from fOffsetCache mutable sbn::pot::BPMOffsetCache_t fOffsetCache; ///< last valid BPM offsets seen + bool fImproveFOM; ///< whether to run sbn::improveBNBqualityFOMs() at the end of the subrun + sbn::BNBFOMFillConfig fFOMFillConfig; ///< its configuration /// Returns the information of the trigger in the current event. sbn::pot::TriggerInfo_t extractTriggerInfo(art::Event const& e) const; @@ -244,6 +262,7 @@ sbn::ICARUSBNBRetriever::ICARUSBNBRetriever(Parameters const& params) bfp_mwr( ifbeam_handle->getBeamFolder(params().MultiWireBundle(), params().URL(), params().MWR_TimeWindow())), fMWRMaxTimeDiff(params().MWRMaxTimeDiff()), fReuseLastBPMOffsets(params().ReuseLastBPMOffsets()), + fImproveFOM(params().ImproveFOM()), fTriggerDatabaseFile(params().TriggerDatabaseFile()) { @@ -267,6 +286,9 @@ sbn::ICARUSBNBRetriever::ICARUSBNBRetriever(Parameters const& params) //bfp_mwr->setValidWindow(86400); bfp_mwr->setValidWindow(3605); produces< std::vector< sbn::BNBSpillInfo >, art::InSubRun >(); + produces< std::vector< unsigned int >, art::InSubRun >("fomStatus"); + fFOMFillConfig.nbMinFOM = params().NeighbourMinFOM(); + fFOMFillConfig.burstFill = params().FillBPMDropouts(); TotalBeamSpills = 0; cet::search_path sp("FW_SEARCH_PATH"); @@ -476,6 +498,13 @@ mf::LogDebug("ICARUSBNBRetriever")<< "Total number of Selected Spills : " << fOu } } + // Beam-quality status of each spill, and the FOM of spills with missing + // inputs from their neighbours (empty multiwire, BPM gaps, 875 drop-outs). + auto status = std::make_unique< std::vector >(); + if (fImproveFOM) *status = sbn::improveBNBqualityFOMs(fOutbeamInfos, fFOMFillConfig); + else for (sbn::BNBSpillInfo const& info: fOutbeamInfos) status->push_back(sbn::getBNBBeamState(info).status); + sr.put(std::move(status), "fomStatus", art::subRunFragment()); + auto p = std::make_unique< std::vector< sbn::BNBSpillInfo > >(); std::swap(*p, fOutbeamInfos); diff --git a/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc index 8830f9cb0..d7bb434ac 100644 --- a/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/SBNDBNBRetriever/SBNDBNBRetriever_module.cc @@ -10,6 +10,7 @@ #include "sbnobj/Common/POTAccounting/BNBSpillInfo.h" #include "sbncode/BeamSpillInfoRetriever/POTTools.h" #include "sbncode/BeamSpillInfoRetriever/getFOM.h" +#include "sbncode/BeamSpillInfoRetriever/BNBFOMFill.h" namespace sbn { class SBNDBNBRetriever; @@ -46,6 +47,8 @@ class sbn::SBNDBNBRetriever : public art::EDProducer { double fMWRMaxTimeDiff; ///< largest accepted multiwire-spill time difference [s] bool fReuseLastBPMOffsets; ///< whether to fill missing BPM offsets from fOffsetCache mutable sbn::pot::BPMOffsetCache_t fOffsetCache; ///< last valid BPM offsets seen + bool fImproveFOM; ///< whether to run sbn::improveBNBqualityFOMs() at the end of the subrun + sbn::BNBFOMFillConfig fFOMFillConfig; ///< its configuration sbn::pot::TriggerInfo_t extractTriggerInfo(art::Event const& e) const; int matchMultiWireData( @@ -60,6 +63,10 @@ class sbn::SBNDBNBRetriever : public art::EDProducer { sbn::SBNDBNBRetriever::SBNDBNBRetriever(fhicl::ParameterSet const & params) : EDProducer{params} { produces< std::vector< sbn::BNBSpillInfo >, art::InSubRun >(); + produces< std::vector< unsigned int >, art::InSubRun >("fomStatus"); + fImproveFOM = params.get("ImproveFOM", true); + fFOMFillConfig.nbMinFOM = params.get("NeighbourMinFOM", -1.); + fFOMFillConfig.burstFill = params.get("FillBPMDropouts", false); fTimePad = params.get("TimePadding"); fBESOffset = params.get("BESOffset"); fDeviceUsedForTiming = params.get("DeviceUsedForTiming"); @@ -230,6 +237,13 @@ void sbn::SBNDBNBRetriever::endSubRun(art::SubRun& sr) } } + // Beam-quality status of each spill, and the FOM of spills with missing + // inputs from their neighbours (empty multiwire, BPM gaps, 875 drop-outs). + auto status = std::make_unique< std::vector >(); + if (fImproveFOM) *status = sbn::improveBNBqualityFOMs(fOutbeamInfos, fFOMFillConfig); + else for (sbn::BNBSpillInfo const& info: fOutbeamInfos) status->push_back(sbn::getBNBBeamState(info).status); + sr.put(std::move(status), "fomStatus", art::subRunFragment()); + auto p = std::make_unique< std::vector< sbn::BNBSpillInfo > >(); std::swap(*p, fOutbeamInfos); diff --git a/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc index d0e88bc25..5279b8eed 100644 --- a/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc +++ b/sbncode/BeamSpillInfoRetriever/SBNDBNBZEROBIASRetriever/SBNDBNBZEROBIASRetriever_module.cc @@ -10,6 +10,7 @@ #include "sbnobj/Common/POTAccounting/BNBSpillInfo.h" #include "sbncode/BeamSpillInfoRetriever/POTTools.h" #include "sbncode/BeamSpillInfoRetriever/getFOM.h" +#include "sbncode/BeamSpillInfoRetriever/BNBFOMFill.h" namespace sbn { class SBNDBNBZEROBIASRetriever; @@ -46,6 +47,8 @@ class sbn::SBNDBNBZEROBIASRetriever : public art::EDProducer { mutable sbn::pot::BPMOffsetCache_t fOffsetCache; ///< last valid BPM offsets seen std::vector< sbn::BNBSpillInfo > fOutbeamInfos; std::vector< sbn::BNBSpillInfo > fOutbeamInfosTotal; + bool fImproveFOM; ///< whether to run sbn::improveBNBqualityFOMs() on the subrun spills + sbn::BNBFOMFillConfig fFOMFillConfig; ///< its configuration sbn::pot::TriggerInfo_t extractTriggerInfo(art::Event const& e) const; void matchMultiWireData( @@ -79,6 +82,10 @@ sbn::SBNDBNBZEROBIASRetriever::SBNDBNBZEROBIASRetriever(fhicl::ParameterSet cons TotalBeamSpills = 0; produces< std::vector< sbn::BNBSpillInfo >, art::InEvent >(); produces< std::vector< sbn::BNBSpillInfo >, art::InSubRun >(); + produces< std::vector< unsigned int >, art::InSubRun >("fomStatus"); + fImproveFOM = params.get("ImproveFOM", true); + fFOMFillConfig.nbMinFOM = params.get("NeighbourMinFOM", -1.); + fFOMFillConfig.burstFill = params.get("FillBPMDropouts", false); } void sbn::SBNDBNBZEROBIASRetriever::produce(art::Event & e) @@ -234,6 +241,13 @@ void sbn::SBNDBNBZEROBIASRetriever::endSubRun(art::SubRun& sr) { mf::LogDebug("SBNDBNBZEROBIASRetriever")<< "Total number of DAQ Spills : " << TotalBeamSpills << std::endl; mf::LogDebug("SBNDBNBZEROBIASRetriever")<< "Total number of Selected Spills : " << fOutbeamInfosTotal.size() << std::endl; + // Beam-quality status of each spill, and the FOM of spills with missing + // inputs from their neighbours (empty multiwire, BPM gaps, 875 drop-outs). + auto status = std::make_unique< std::vector >(); + if (fImproveFOM) *status = sbn::improveBNBqualityFOMs(fOutbeamInfosTotal, fFOMFillConfig); + else for (sbn::BNBSpillInfo const& info: fOutbeamInfosTotal) status->push_back(sbn::getBNBBeamState(info).status); + sr.put(std::move(status), "fomStatus", art::subRunFragment()); + auto p = std::make_unique< std::vector< sbn::BNBSpillInfo > >(); std::swap(*p, fOutbeamInfosTotal); sr.put(std::move(p), art::subRunFragment()); diff --git a/sbncode/BeamSpillInfoRetriever/getFOM.cpp b/sbncode/BeamSpillInfoRetriever/getFOM.cpp index c7bebfd90..a22dcf0dc 100644 --- a/sbncode/BeamSpillInfoRetriever/getFOM.cpp +++ b/sbncode/BeamSpillInfoRetriever/getFOM.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include "TH1D.h" #include "TFitResult.h" #include @@ -43,113 +44,214 @@ namespace sbn { } // local namespace - std::tuple getBNBqualityFOM(BNBSpillInfo const& spill) - { - //Z Position of the monitors in m - double const vp873_zpos= 191.153656; - double const hp875_zpos= 202.116104; - double const vp875_zpos= 202.3193205; - double const hptg1_zpos= 204.833267; - double const hptg2_zpos= 205.240662; - double const vptg2_zpos= 205.036835; - double const target_center_zpos= 206.870895; - double const p875x[]={0.431857, 0.158077, 0.00303551}; - double const p875y[]={0.279128, 0.337048, 0}; - double const p876x[]={0.166172, 0.30999, -0.00630299}; - double const p876y[]={0.13425, 0.580862, 0}; - - // ---- intensity: TOR860, falling back to TOR875. - // A missing toroid is stored as -999e12; a spill with no (or negative) - // intensity has no meaningful FOM (the optics model needs ppp > 0). - double tor = -1.; - if (isValid(spill.TOR860) && spill.TOR860 > 0.) tor = spill.TOR860; - else if (isValid(spill.TOR875) && spill.TOR875 > 0.) tor = spill.TOR875; - else return {-1, -1, -1}; - - // ---- horizontal: HP875 and HPTG1 (HPTG2 as fallback), offsets subtracted. - // Missing devices used to be fed into the extrapolation as -999 (the - // previous `.empty()` checks could never trigger), placing the beam ~1 m - // off target and losing the spill. - bool const okHP875 = isValid(spill.HP875) && isValid(spill.HP875Offset); - bool const okHPTG1 = isValid(spill.HPTG1) && isValid(spill.HPTG1Offset); - bool const okHPTG2 = isValid(spill.HPTG2) && isValid(spill.HPTG2Offset); - if (!okHP875 || (!okHPTG1 && !okHPTG2)) return {2, 2, 2}; - double const delta_hp875 = spill.HP875 - spill.HP875Offset; - auto const [ horang, horpos ] = okHPTG1 - ? extrapolate(delta_hp875, hp875_zpos, spill.HPTG1 - spill.HPTG1Offset, hptg1_zpos, target_center_zpos) - : extrapolate(delta_hp875, hp875_zpos, spill.HPTG2 - spill.HPTG2Offset, hptg2_zpos, target_center_zpos); - - // ---- vertical: VP875 and VP873 (VPTG2 as fallback), offsets subtracted. - bool const okVP875 = isValid(spill.VP875) && isValid(spill.VP875Offset); - bool const okVP873 = isValid(spill.VP873) && isValid(spill.VP873Offset); - bool const okVPTG2 = isValid(spill.VPTG2) && isValid(spill.VPTG2Offset); - if (!okVP875 || (!okVP873 && !okVPTG2)) return {3, 3, 3}; - double const delta_vp875 = spill.VP875 - spill.VP875Offset; - auto const [ verang, verpos ] = okVP873 - ? extrapolate(delta_vp875, vp875_zpos, spill.VP873 - spill.VP873Offset, vp873_zpos, target_center_zpos) - : extrapolate(delta_vp875, vp875_zpos, spill.VPTG2 - spill.VPTG2Offset, vptg2_zpos, target_center_zpos); - - const double smallSigmaX =0.5, largeSigmaX = 10, smallSigmaY = 0.3, largeSigmaY =10, maxChi2X = 20, maxChi2Y = 20; - auto inWindow = [&](double sx, double sy) - { return sx>smallSigmaX && sxsmallSigmaY && sy 0`, - // allowed reading past the end of a short vector) and both fits succeed. - // Preference: target multiwire (no transformation), then M876, then M875. - struct MWDevice_t { - std::vector const* data; - double const* px; double const* py; - }; - MWDevice_t const mwDevices[] = { - { &spill.MMBTBB, nullptr, nullptr }, - { &spill.M876BB, p876x, p876y }, - { &spill.M875BB, p875x, p875y }, - }; - double tgtsx = MissingValue, tgtsy = MissingValue; - bool goodFit = false; - for (auto const& dev: mwDevices) { - if (dev.data->size() < 2*NWires) continue; - std::vector const mw(dev.data->begin(), dev.data->begin() + 2*NWires); - double xx, yy, sx, sy, chi2x, chi2y; - bool const fitOK = processBNBprofile(&mw[0], xx, sx, chi2x) - & processBNBprofile(&mw[NWires], yy, sy, chi2y); - if (!fitOK) continue; - if (dev.px) { - sx = dev.px[0] + dev.px[1]*sx + dev.px[2]*sx*sx; - sy = dev.py[0] + dev.py[1]*sy + dev.py[2]*sy*sy; + namespace { + + // Z positions of the monitors [m] + constexpr double vp873_zpos = 191.153656; + constexpr double hp875_zpos = 202.116104; + constexpr double vp875_zpos = 202.3193205; + constexpr double hptg1_zpos = 204.833267; + constexpr double hptg2_zpos = 205.240662; + constexpr double vptg2_zpos = 205.036835; + constexpr double target_center_zpos = 206.870895; + + // multiwire sigma -> sigma at the target (quadratic polynomials) + constexpr double p875x[] = {0.431857, 0.158077, 0.00303551}; + constexpr double p875y[] = {0.279128, 0.337048, 0}; + constexpr double p876x[] = {0.166172, 0.30999, -0.00630299}; + constexpr double p876y[] = {0.13425, 0.580862, 0}; + + constexpr double smallSigmaX = 0.5, largeSigmaX = 10, smallSigmaY = 0.3, largeSigmaY = 10; + constexpr double maxChi2X = 20, maxChi2Y = 20; + bool inWindow(double sx, double sy) + { return sx > smallSigmaX && sx < largeSigmaX && sy > smallSigmaY && sy < largeSigmaY; } + + /// A database (online-fitted) multiwire sigma above this [mm] means the + /// chamber saw no beam: the online fit went to noise (beam: 1-2.5 mm, empty + /// profile: ~7 mm). Such a width is not used. + constexpr double EmptyMWSigma = 4.0; + + double poly(double const* p, double s) { return p[0] + p[1]*s + p[2]*s*s; } + + /** + * Autotune calibration times [s, UTC] of the two horizontal target BPMs + * (from IFBeam, Z. Pavlovic; same table as sbnana's getBNBFoM.cxx). + * Autotune steers the beam on HP875 and on the target BPM calibrated most + * recently, so that one reads the beam at its setpoint; the offset of the + * other one may be stale (HPTG1 before 2023-03-03: +1.5 mm, i.e. the beam + * projected 2.6 mm off-centre). New autotune calibrations must be added here. + */ + constexpr unsigned long HPTG1CalibTimes[] = { 1420092000, 1576014670, 1606510357, 1677887110, 1711030691 }; + constexpr unsigned long HPTG2CalibTimes[] = { 1420092000, 1574190444, 1576014670, 1588603117, 1605706657, + 1606510357, 1608593857, 1668802510 }; + + template + unsigned long lastCalibration(unsigned long const (×)[N], unsigned long t) { + unsigned long last = 0; + for (unsigned long c: times) if (c <= t) last = c; + return last; + } + + /// Fills position, angle and the related status bits of `state`. + void fillGeometry(BNBSpillInfo const& spill, BNBBeamState& state) { + using namespace fomstatus; + // ---- horizontal: HP875 and the primary target BPM, offsets subtracted. + // The other target BPM is not used as a fallback: its offset is not the + // one autotune steers on, so the position would be biased (mm level). + bool const useHPTG2 = hptg2IsPrimary(spill.spill_time_s); + if (useHPTG2) state.status |= HPTG2Primary; + bool const okHP875 = isValid(spill.HP875) && isValid(spill.HP875Offset); + bool const okHPTG1 = isValid(spill.HPTG1) && isValid(spill.HPTG1Offset); + bool const okHPTG2 = isValid(spill.HPTG2) && isValid(spill.HPTG2Offset); + bool const okPrimary = useHPTG2? okHPTG2: okHPTG1; + if (!okHP875 || !okPrimary) { + state.status |= NoHBPM; + if (okHP875 && (useHPTG2? okHPTG1: okHPTG2)) state.status |= FallbackRemoved; + } + else { + double const delta_hp875 = spill.HP875 - spill.HP875Offset; + std::tie(state.hang, state.hpos) = useHPTG2 + ? extrapolate(delta_hp875, hp875_zpos, spill.HPTG2 - spill.HPTG2Offset, hptg2_zpos, target_center_zpos) + : extrapolate(delta_hp875, hp875_zpos, spill.HPTG1 - spill.HPTG1Offset, hptg1_zpos, target_center_zpos); } - if (inWindow(sx, sy) && chi2x < maxChi2X && chi2y < maxChi2Y) { - tgtsx = sx; tgtsy = sy; - goodFit = true; - break; + + // ---- vertical: VP875 and VP873, offsets subtracted (VPTG2 is not used + // as a fallback for the same reason). + bool const okVP875 = isValid(spill.VP875) && isValid(spill.VP875Offset); + bool const okVP873 = isValid(spill.VP873) && isValid(spill.VP873Offset); + bool const okVPTG2 = isValid(spill.VPTG2) && isValid(spill.VPTG2Offset); + if (!okVP875 || !okVP873) { + state.status |= NoVBPM; + if (okVP875 && okVPTG2) state.status |= FallbackRemoved; + } + else { + double const delta_vp875 = spill.VP875 - spill.VP875Offset; + std::tie(state.vang, state.vpos) = + extrapolate(delta_vp875, vp875_zpos, spill.VP873 - spill.VP873Offset, vp873_zpos, target_center_zpos); } } - double const fom = goodFit - ? 1-pow(10,sbn::calcFOM(horpos,horang,verpos,verang,tor,tgtsx,tgtsy)) - : MissingValue; - - // ---- "pre-fit" FOM with the widths fitted online (database M876/M875). - // There is no chi2 for these: the old code applied the chi2 of whichever - // multiwire profile it fitted last (uninitialised if none was fitted). - double prefitfom = MissingValue; - struct DBWidth_t { double hs, vs; double const* px; double const* py; }; - DBWidth_t const dbWidths[] = { - { spill.M876HS, spill.M876VS, p876x, p876y }, - { spill.M875HS, spill.M875VS, p875x, p875y }, - }; - for (auto const& w: dbWidths) { - if (!isValid(w.hs) || !isValid(w.vs)) continue; - double const sx = w.px[0] + w.px[1]*w.hs + w.px[2]*w.hs*w.hs; - double const sy = w.py[0] + w.py[1]*w.vs + w.py[2]*w.vs*w.vs; - if (!inWindow(sx, sy)) continue; - prefitfom = 1-pow(10,sbn::calcFOM(horpos,horang,verpos,verang,tor,sx,sy)); - break; + + /// Width from the multiwire profiles fitted here (M876, then M875). + /// The target multiwire (MMBTBB) is not used: its profiles are noise. + bool fittedWidth(BNBSpillInfo const& spill, double& tgtsx, double& tgtsy, int& source) { + struct MWDevice_t { std::vector const* data; double const* px; double const* py; int source; }; + MWDevice_t const mwDevices[] = { + { &spill.M876BB, p876x, p876y, BNBBeamState::FitM876 }, + { &spill.M875BB, p875x, p875y, BNBBeamState::FitM875 }, + }; + for (auto const& dev: mwDevices) { + // each profile is 48 horizontal wires followed by 48 vertical ones + if (dev.data->size() < 2*NWires) continue; + std::vector const mw(dev.data->begin(), dev.data->begin() + 2*NWires); + double xx, yy, sx, sy, chi2x, chi2y; + bool const fitOK = processBNBprofile(&mw[0], xx, sx, chi2x) + & processBNBprofile(&mw[NWires], yy, sy, chi2y); + if (!fitOK) continue; + sx = poly(dev.px, sx); + sy = poly(dev.py, sy); + if (inWindow(sx, sy) && chi2x < maxChi2X && chi2y < maxChi2Y) { + tgtsx = sx; tgtsy = sy; source = dev.source; + return true; + } + } + return false; + } + + /// Width from the online ("pre-fit", database) multiwire widths. + /// Sets `empty` if the device chosen saw no beam (no width returned then). + bool databaseWidth(BNBSpillInfo const& spill, double& tgtsx, double& tgtsy, int& source, bool& empty) { + struct DBWidth_t { double hs, vs; double const* px; double const* py; int source; }; + DBWidth_t const dbWidths[] = { + { spill.M876HS, spill.M876VS, p876x, p876y, BNBBeamState::DatabaseM876 }, + { spill.M875HS, spill.M875VS, p875x, p875y, BNBBeamState::DatabaseM875 }, + }; + empty = false; + for (auto const& w: dbWidths) { + if (!isValid(w.hs) || !isValid(w.vs)) continue; + double const sx = poly(w.px, w.hs); + double const sy = poly(w.py, w.vs); + if (!inWindow(sx, sy)) continue; + if (w.hs > EmptyMWSigma || w.vs > EmptyMWSigma) { empty = true; return false; } + tgtsx = sx; tgtsy = sy; source = w.source; + return true; + } + return false; } - // ---- FOM with the nominal beam width (scale factors 1) - double const noMWfom = 1-pow(10,sbn::calcFOM(horpos,horang,verpos,verang,tor)); + } // local namespace + + + bool hptg2IsPrimary(unsigned long spill_time_s) { + return lastCalibration(HPTG2CalibTimes, spill_time_s) >= lastCalibration(HPTG1CalibTimes, spill_time_s); + } + + + double beamIntensity(BNBSpillInfo const& spill, bool* tor875used) { + // TOR860 sometimes reads ~1e-3 of the beam TOR875 sees (TOR860 drop-out): + // then TOR875 is the intensity. A missing toroid is -999. + bool const ok860 = isValid(spill.TOR860) && spill.TOR860 > 0.; + bool const ok875 = isValid(spill.TOR875) && spill.TOR875 > 0.; + bool const dropout = ok860 && ok875 && spill.TOR860 < TOR860DropoutFraction * spill.TOR875; + if (tor875used) *tor875used = dropout; + if (ok860 && !dropout) return spill.TOR860; + if (ok875) return spill.TOR875; + return -1.; + } + + + BNBBeamState getBNBBeamState(BNBSpillInfo const& spill) { + using namespace fomstatus; + BNBBeamState state; + bool tor875used = false; + state.tor = beamIntensity(spill, &tor875used); + if (tor875used) state.status |= TOR875Used; + if (state.tor <= 0.) state.status |= NoTOR; + fillGeometry(spill, state); + + // width: fitted profile, else database width, else nominal + bool empty = false; + if (!fittedWidth(spill, state.sx, state.sy, state.widthSource) + && !databaseWidth(spill, state.sx, state.sy, state.widthSource, empty)) + { + state.sx = state.sy = MissingValue; + state.widthSource = BNBBeamState::Nominal; + state.status |= NoMWWidth; + if (empty) state.status |= MWEmpty; + } + return state; + } + + + double computeFOM(BNBBeamState const& state, bool nominalWidth) { + if (!state.hasFOM()) return MissingValue; + bool const measured = !nominalWidth && state.widthSource != BNBBeamState::Nominal; + return measured + ? 1-pow(10, sbn::calcFOM(state.hpos, state.hang, state.vpos, state.vang, state.tor, state.sx, state.sy)) + : 1-pow(10, sbn::calcFOM(state.hpos, state.hang, state.vpos, state.vang, state.tor)); + } + + + std::tuple getBNBqualityFOM(BNBSpillInfo const& spill) + { + BNBBeamState const state = getBNBBeamState(spill); + if (state.status & fomstatus::NoTOR) return {-1, -1, -1}; + if (state.status & fomstatus::NoHBPM) return {2, 2, 2}; + if (state.status & fomstatus::NoVBPM) return {3, 3, 3}; + + // FOM with the width fitted from the multiwire profiles + double fom = MissingValue, prefitfom = MissingValue; + double sx = MissingValue, sy = MissingValue; + int source = BNBBeamState::Nominal; + if (fittedWidth(spill, sx, sy, source)) + fom = 1-pow(10, sbn::calcFOM(state.hpos, state.hang, state.vpos, state.vang, state.tor, sx, sy)); + + // "pre-fit" FOM with the widths fitted online (database M876/M875) + bool empty = false; + if (databaseWidth(spill, sx, sy, source, empty)) + prefitfom = 1-pow(10, sbn::calcFOM(state.hpos, state.hang, state.vpos, state.vang, state.tor, sx, sy)); + + // FOM with the nominal beam width (scale factors 1) + double const noMWfom = computeFOM(state, true); return {fom, prefitfom, noMWfom}; } diff --git a/sbncode/BeamSpillInfoRetriever/getFOM.h b/sbncode/BeamSpillInfoRetriever/getFOM.h index a3cba4ac9..e83cf02c7 100644 --- a/sbncode/BeamSpillInfoRetriever/getFOM.h +++ b/sbncode/BeamSpillInfoRetriever/getFOM.h @@ -15,6 +15,55 @@ namespace sbn { + /// Bits of the beam-quality status word (same values as the analysis-level + /// `fom_status` of bnb_fom.py / fom_fixups.py). + namespace fomstatus { + constexpr unsigned int NoTOR = 1; ///< no valid, positive TOR860 or TOR875 (no beam) + constexpr unsigned int NoHBPM = 2; ///< horizontal position/angle cannot be formed + constexpr unsigned int NoVBPM = 4; ///< vertical position/angle cannot be formed + constexpr unsigned int NoMWWidth = 64; ///< no usable multiwire width: nominal width + constexpr unsigned int NeighborFilled = 128; ///< position/angle from the adjacent spills + constexpr unsigned int TOR875Used = 256; ///< TOR860 drop-out: intensity from TOR875 + constexpr unsigned int MWEmpty = 512; ///< database multiwire width is an empty chamber + constexpr unsigned int WidthFromNeighbors = 1024; ///< width from the adjacent spills + constexpr unsigned int FallbackRemoved = 2048; ///< only the non-primary BPM read: not used + constexpr unsigned int BurstFill = 16384; ///< 875-station drop-out filled from both sides + constexpr unsigned int HPTG2Primary = 32768; ///< horizontal projection through HPTG2 + } + + /// TOR860 below this fraction of TOR875 is a TOR860 drop-out. + constexpr double TOR860DropoutFraction = 0.5; + + /// Beam at the centre of the target, as used by the figure of merit. + struct BNBBeamState { + enum WidthSource_t: int { FitM876 = 1, FitM875 = 2, DatabaseM876 = 3, DatabaseM875 = 4, Nominal = 5 }; + double tor = -1.; ///< intensity used [protons] + double hpos = -999.; ///< horizontal position [mm] + double hang = -999.; ///< horizontal angle [mrad] + double vpos = -999.; ///< vertical position [mm] + double vang = -999.; ///< vertical angle [mrad] + double sx = -999.; ///< horizontal width at the target [mm] (-999: nominal) + double sy = -999.; ///< vertical width at the target [mm] (-999: nominal) + int widthSource = Nominal; + unsigned int status = 0; ///< bits from `sbn::fomstatus` + bool hasFOM() const + { return (status & (fomstatus::NoTOR | fomstatus::NoHBPM | fomstatus::NoVBPM)) == 0; } + }; + + /// Whether HPTG2 (instead of HPTG1) is the primary horizontal target BPM at + /// this time: the one autotune was calibrated on most recently. + bool hptg2IsPrimary(unsigned long spill_time_s); + + /// Intensity for the FOM [protons]: TOR860, or TOR875 if TOR860 is missing or + /// dropped out (`*tor875used` is set in the latter case); -1 if none. + double beamIntensity(BNBSpillInfo const& spill, bool* tor875used = nullptr); + + /// Position, angle, width and status of the beam at the target for one spill. + BNBBeamState getBNBBeamState(BNBSpillInfo const& spill); + + /// FOM of a beam state (-999 if it has none); `nominalWidth` ignores its width. + double computeFOM(BNBBeamState const& state, bool nominalWidth = false); + /** * @brief Returns a Figure of Merit on BNB beam quality. * @@ -30,6 +79,12 @@ namespace sbn * * `3`: vertical position/angle cannot be formed (missing BPM or BPM offset) * The first two values are additionally `-999` when no usable width is found. * A device value of `-999` is treated as missing. + * + * Horizontal: HP875 and the primary target BPM (`hptg2IsPrimary()`); + * vertical: VP875 and VP873. The other target BPMs are not used as fallbacks. + * Widths: M876 then M875 profile fits (not the target multiwire), then the + * database widths unless the chamber was empty (sigma > 4 mm). + * Intensity: `beamIntensity()`. */ std::tuple getBNBqualityFOM(BNBSpillInfo const& spill); diff --git a/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl b/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl index 30d39e36c..db1b2dbfa 100644 --- a/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl +++ b/sbncode/BeamSpillInfoRetriever/job/icarusbnbspillinfo.fcl @@ -18,5 +18,8 @@ icarusbnbspillinfo: { TriggerDatabaseFile: "triggerDatabase/icarus_triggers.db" MWRMaxTimeDiff: 0.0333 #unit seconds, multiwire readings further than this from the spill are not used (<= 0 disables) ReuseLastBPMOffsets: true #if a BPM offset query fails, use the last valid offset seen in this job + ImproveFOM: true #recover the FOM of spills with missing inputs from their neighbours (BNBFOMFill.h) + NeighbourMinFOM: 0.998 #neighbour fill only between spills with FOM > 0.998 (closure test, ICARUS Run 2) + FillBPMDropouts: true #fill HP875+VP875 drop-outs from the target BPMs + measured spills on both sides } END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl b/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl index 524586264..246af4716 100644 --- a/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl +++ b/sbncode/BeamSpillInfoRetriever/job/sbndbnbdefaults.fcl @@ -16,6 +16,9 @@ sbndbnbspillinfo: { DeviceUsedForTiming: "E:TOR860" MWRMaxTimeDiff: 0.0333 #unit seconds, multiwire readings further than this from the spill are not used (<= 0 disables) ReuseLastBPMOffsets: true #if a BPM offset query fails, use the last valid offset seen in this job + ImproveFOM: true #recover the FOM of spills with missing inputs from their neighbours (BNBFOMFill.h) + NeighbourMinFOM: -1 #no FOM requirement on the neighbours (SBND closure test: 43 bad spills in 909k) + FillBPMDropouts: false #875 drop-outs in SBND are single spills, covered by the neighbour fill } END_PROLOG diff --git a/sbncode/CAFMaker/FillExposure.cxx b/sbncode/CAFMaker/FillExposure.cxx index 12decbfa8..5d656a489 100644 --- a/sbncode/CAFMaker/FillExposure.cxx +++ b/sbncode/CAFMaker/FillExposure.cxx @@ -26,9 +26,12 @@ namespace caf } else if((NoMultiWireFOM >= 0.0) && (NoMultiWireFOM <= 1.0)) { - finalFOM = 100.+NoMultiWireFOM; + // nominal-width FOM: used as is (it used to be stored as 100 + FOM, + // which made the standard 0.98 < FOM <= 1 cut reject these spills; + // the beam-quality status, product "fomStatus" of the retriever, has + // bit 64 set for them) + finalFOM = NoMultiWireFOM; std::cout << "makeSRBNBInfo: Chose assumed width of 1.0" << std::endl; - std::cout << "makeSRBNBInfo: Note, that this gets prefixed with 100 so it doesn't get automatically included in standard FOM cuts" << std::endl; } caf::SRBNBInfo single_store;