diff --git a/CMakeLists.txt b/CMakeLists.txt index 6e268929e..0c4c211d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -15,7 +15,7 @@ cmake_minimum_required(VERSION 3.19 FATAL_ERROR) -project(sbncode VERSION 09.32.01 LANGUAGES CXX) +project(sbncode VERSION 09.32.01.01 LANGUAGES CXX) message(STATUS "\n\n ========================== ${PROJECT_NAME} ==========================") @@ -64,14 +64,19 @@ find_ups_product(larreco) find_ups_product(larpandora) find_ups_product(larsoft v09_00_00) find_ups_product(artdaq_core v3_06_01) +find_ups_product(ifbeam) +find_ups_product(ifbeam_art) +find_ups_product(libwda) find_ups_product(sbnobj v09_10_00) find_ups_product(sbndata v01_00) find_ups_product(ifdhc) find_ups_product(ifdh_art) find_ups_product(log4cpp) +find_ups_boost() find_ups_root() #find_ups_product(sbndcode v06_67_00) #find_ups_product(uboonecode v06_67_00) +find_ups_product( sbndaq_artdaq_core ) # macros for dictionary and simple_plugin include(ArtDictionary) diff --git a/sbncode/BeamSpillInfoRetriever/BNBRetriever/BNBRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/BNBRetriever/BNBRetriever_module.cc new file mode 100644 index 000000000..9e954f1b2 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/BNBRetriever/BNBRetriever_module.cc @@ -0,0 +1,600 @@ +/** ******************************************************************** + * @file BNBRetriever_module.cc + * @date Wed April 9 2021 + * @author J. Zennamo (FNAL) + * + * Based heavily on code by Z. Pavlovic written for MicroBooNE + * Based heavily on code by NOvA collaboration (Thanks NOvA!): + * https://cdcvs.fnal.gov/redmine/projects/novaart/repository/entry/trunk/IFDBSpillInfo/BNBInfo_module.cc + */ + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "larcorealg/CoreUtils/counter.h" + +#include "artdaq-core/Data/Fragment.hh" +#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerUDPFragment.hh" + +#include "sbnobj/Common/POTAccounting/BNBSpillInfo.h" + +#include "IFBeam_service.h" +#include "ifbeam_c.h" +#include "MWRData.h" + +#include +#include +#include + +namespace sbn { + class BNBRetriever; +} + +class sbn::BNBRetriever : public art::EDProducer { +public: + + struct Config { + + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Atom TimePadding { + Name{ "TimePadding" }, + Comment{ "extension to the time window considered when collecting spills [seconds]" }, + 0.0333 // default + }; + + fhicl::Atom RawDataLabel { + Name{ "raw_data_label" }, + Comment{ "art data product instance name for trigger information (product label is 'daq')" } + }; + + fhicl::Atom DeviceUsedForTiming { + Name{ "DeviceUsedForTiming" }, + Comment{ "name in the IFBeam database of the device used to extract spill times" } + }; + + fhicl::Atom URL { + Name{ "URL" }, + Comment{ "IFBeam database access URL" } + }; + + fhicl::Atom Bundle { + Name{ "Bundle" }, + Comment{ "" } // explain what this is and which database/table it's looking for + }; + + fhicl::Atom TimeWindow { + Name{ "TimeWindow" }, + Comment{ "" } // explain what this is, what's for and its unit + }; + + fhicl::Atom MultiWireBundle { + Name{ "MultiWireBundle" }, + Comment{ "" } // explain what this is and which database/table it's looking for + }; + + fhicl::Atom MWR_TimeWindow { + Name{ "MWR_TimeWindow" }, + Comment{ "" } // explain what this is, what's for and its unit + }; + + }; // Config + + using Parameters = art::EDProducer::Table; + + + explicit BNBRetriever(Parameters const& params); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + BNBRetriever(BNBRetriever const&) = delete; + BNBRetriever(BNBRetriever&&) = delete; + BNBRetriever& operator=(BNBRetriever const&) = delete; + BNBRetriever& operator=(BNBRetriever&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + void beginSubRun(art::SubRun& sr) override; + void endSubRun(art::SubRun& sr) override; + +private: + // input labels + std::vector< sbn::BNBSpillInfo > fOutbeamInfos; + double fTimePad; + std::string fURL; + MWRData mwrdata; + std::string raw_data_label; + std::string fDeviceUsedForTiming; + unsigned int TotalBeamSpills; + // + art::ServiceHandle ifbeam_handle; + std::unique_ptr bfp; + std::unique_ptr bfp_mwr; + + struct TriggerInfo_t { + int gate_type = 0; ///< Source of the spill: `1`: BNB, `2`: NuMI + double t_current_event = 0; + double t_previous_event = 0; + unsigned int number_of_gates_since_previous_event = 0; // FIXME needs to be integral type + }; + + struct MWRdata_t { + std::vector< std::vector > MWR_times; + std::vector< std::vector< std::vector< int > > > unpacked_MWR; + }; + + + static constexpr double MWRtoroidDelay = -0.035; ///< the same time point is measured _t_ by MWR and _t + MWRtoroidDelay`_ by the toroid [ms] + + /// Returns the information of the trigger in the current event. + TriggerInfo_t extractTriggerInfo(art::Event const& e) const; + + /** + * @brief Determines spill times and extracts data based on multiwire devices. + * @param triggerInfo information from the trigger of this event + * @return times and unpacked data, per device (`"M875BB"`, `"M876BB"`, `"MMBTBB"`) + */ + MWRdata_t extractSpillTimes(TriggerInfo_t const& triggerInfo) const; + + /** + * @brief Matches spill times with multiwire chamber data from the database. + * @param triggerInfo information from the trigger of this event + * @param MWRdata data from multiwire chambers + * @param isFirstEventInRun whether we are processing the first event of the run + * @param[out] beamInfos container to _add_ spill information records to + * @return count of matched spills + */ + int matchMultiWireData( + TriggerInfo_t const& triggerInfo, + MWRdata_t const& MWRdata, bool isFirstEventInRun, + std::vector< sbn::BNBSpillInfo >& beamInfos + ) const; + + /** + * @brief Assembles and returns a spill information record. + * @param time time of the spill + * @param MWRdata all extracted data from multiwire chambers + * @param matched_MWR data from multiwire chambers matched with the time + * @return a `sbn::BNBSpillInfo` object with information on the spill at `time` + */ + sbn::BNBSpillInfo makeBNBSpillInfo + (double time, MWRdata_t const& MWRdata, std::vector const& matched_MWR) const; + +}; + +sbn::BNBRetriever::BNBRetriever(Parameters const& params) + : EDProducer{params}, + fTimePad(params().TimePadding()), + raw_data_label(params().RawDataLabel()), + fDeviceUsedForTiming(params().DeviceUsedForTiming()), + bfp( ifbeam_handle->getBeamFolder(params().Bundle(), params().URL(), params().TimeWindow())), + bfp_mwr( ifbeam_handle->getBeamFolder(params().MultiWireBundle(), params().URL(), params().MWR_TimeWindow())) +{ + + // Check fTimePad is positive + if (fTimePad < 0) { + throw art::Exception(art::errors::Configuration) + << "Parameter `TimePadding` must be non-negative (" << fTimePad << " was specified).\n"; + }//End Time check + + // how close in time does the spill time have to be from the DAQ time (in seconds). + // If these are too large then it fails to capture the device + // If these are too small then the time jitter in devices means we miss good data + // + // These values should likely not be changed unless authors of the IFBeam API are consulted + // + bfp->set_epsilon(0.02); //20 ms, this was tuned by hand and compared to IFBeamDB times + bfp_mwr->set_epsilon(0.5); + + //bfp_mwr->setValidWindow(86400); + bfp_mwr->setValidWindow(3605); + produces< std::vector< sbn::BNBSpillInfo >, art::InSubRun >(); + TotalBeamSpills = 0; +} + + +void sbn::BNBRetriever::produce(art::Event& e) +{ + + TriggerInfo_t const triggerInfo = extractTriggerInfo(e); + + //We only want to process BNB gates, i.e. type 1 + if(triggerInfo.gate_type != 1) return; + // Keep track of the number of beam gates the DAQ thinks + // are in this job + TotalBeamSpills += triggerInfo.number_of_gates_since_previous_event; + + + MWRdata_t const MWRdata = extractSpillTimes(triggerInfo); + + + int const spill_count = matchMultiWireData(triggerInfo, MWRdata, e.event() == 1, fOutbeamInfos); + + + if(spill_count > int(triggerInfo.number_of_gates_since_previous_event)) + mf::LogDebug("BNBRetriever")<< "Event Spills : " << spill_count << ", DAQ Spills : " << triggerInfo.number_of_gates_since_previous_event << " \t \t ::: WRONG!"<< std::endl; + else + mf::LogDebug("BNBRetriever")<< "Event Spills : " << spill_count << ", DAQ Spills : " << triggerInfo.number_of_gates_since_previous_event << std::endl; + +}//end iteration over art::Events + + +sbn::BNBRetriever::TriggerInfo_t sbn::BNBRetriever::extractTriggerInfo(art::Event const& e) const { + + //Here we read in the artdaq Fragments and extract three pieces of information: + // 1. The time of the current event, t_current_event + // 2. the time of the previously triggered event, t_previous_event (NOTE: Events are non-sequential!) + // 3. the number of beam spills since the previously triggered event, number_of_gates_since_previous_event + + auto const & raw_data = e.getProduct< std::vector >({ raw_data_label, "ICARUSTriggerUDP" }); + + TriggerInfo_t triggerInfo; + + for(auto raw_datum : raw_data){ + + uint64_t artdaq_ts = raw_datum.timestamp(); + icarus::ICARUSTriggerUDPFragment frag(raw_datum); + std::string data = frag.GetDataString(); + char *buffer = const_cast(data.c_str()); + icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerString(buffer); + triggerInfo.gate_type = datastream_info.gate_type; + triggerInfo.number_of_gates_since_previous_event = frag.getDeltaGatesBNB(); + + triggerInfo.t_current_event = static_cast(artdaq_ts)/(1000000000.0); //check this offset... + if(triggerInfo.gate_type == 1) + triggerInfo.t_previous_event = (static_cast(frag.getLastTimestampBNB()))/(1e9); + else + triggerInfo.t_previous_event = (static_cast(frag.getLastTimestampOther()))/(1000000000.0); + + } + + mf::LogDebug("BNBRetriever") << std::setprecision(19) << "Previous : " << triggerInfo.t_previous_event << ", Current : " << triggerInfo.t_current_event << std::endl; + + return triggerInfo; +} + + +sbn::BNBRetriever::MWRdata_t sbn::BNBRetriever::extractSpillTimes(TriggerInfo_t const& triggerInfo) const { + + // These lines get everything primed within the IFBeamDB + // They seem redundant but they are needed + try{auto cur_vec_temp = bfp->GetNamedVector((triggerInfo.t_previous_event)-fTimePad,"E:THCURR");} catch (WebAPIException &we) {} + try{auto packed_M876BB_temp = bfp_mwr->GetNamedVector((triggerInfo.t_current_event)+fTimePad,"E:M875BB{4440:888}.RAW");} catch (WebAPIException &we) {} + + + //The multiwire chambers provide their + // data in a vector format but we'll have + // to sort through it in std::string format + // to correctly unpack it + std::vector< std::vector< std::vector< int > > > unpacked_MWR; + std::vector< std::vector< double> > MWR_times; + unpacked_MWR.resize(3); + MWR_times.resize(3); + std::string packed_data_str; + + //Create a list of all the MWR devices with their different + // memory buffer increments + // generally in the format: "E:.{Memory Block}" + std::vector vars = bfp_mwr->GetDeviceList(); + mf::LogDebug("BNBRetriever") << " Number of MWR Device Blocks Found : " << vars.size() << std::endl; + // Tracking the time from the IFBeamDB + double time_for_mwr; + + // this is an iterator to track which of the + // three devices we will be working with + int dev = 0; + + // The MWR devices are annoying and have confusing buffer + // what we'll do is sort through all of them first and then + // match them to the closest spills in time + // + + // int t_steps = int(((triggerInfo.t_previous_event - fTimePad) - (triggerInfo.t_current_event + fTimePad))/0.5)+25; + int t_steps = int(((triggerInfo.t_current_event + fTimePad) - (triggerInfo.t_previous_event - fTimePad - 20.))/0.5)+25; + mf::LogDebug("BNBRetriever") << " t_steps " << t_steps << std::endl; + + for(int t = 0; t < t_steps; t++){//Iterate through time increments + for (std::string const& var : vars) {// Iterate through the devices + + //Make sure we have a device + if(var.empty()){ + mf::LogDebug("BNBRetriever") << " NO MWR DEVICES?!" << std::endl; + continue; + } + /// Check the device name and interate the double-vector index + if(var.find("M875BB") != std::string::npos ) dev = 0; + else if(var.find("M876BB") != std::string::npos ) dev = 1; + else if(var.find("MMBTBB") != std::string::npos ) dev = 2; + else{ + mf::LogDebug("BNBRetriever") << " NOT matched to a MWR DEVICES?!" << var << std::endl; + continue;} + + time_for_mwr = 0; + + try{ + //Pull the MWR data for the device + // these data are "packed" + std::vector packed_MWR = bfp_mwr->GetNamedVector((triggerInfo.t_previous_event)-20.-fTimePad+double(0.5*t),var,&time_for_mwr); + + //We'll convert this into a format + // that we can unpack doubles >> strings + // + packed_data_str.clear(); + packed_data_str += std::to_string(int(time_for_mwr)); + packed_data_str.append(","); + packed_data_str.append(var); + packed_data_str.append(",,"); + + /* for(auto const value: packed_MWR){ + packed_data_str += ','; + packed_data_str += std::to_string(int(value)); + }*/ + for(int j = 0; j < int(packed_MWR.size()); j++){ + packed_data_str += std::to_string(int(packed_MWR[j])); + if(j < int(packed_MWR.size())-1) + packed_data_str.append(","); + } + + // Use Zarko's unpacking function to turn this into consumeable data + std::vector MWR_times_temp; + + // There is a 35 ms offset between the toriod and the MWR times + // we'll just remove that here to match to the spill times + std::vector< std::vector< int > > unpacked_MWR_temp = mwrdata.unpackMWR(packed_data_str,MWR_times_temp,MWRtoroidDelay); + + //There are four events that are packed into one MWR IFBeam entry + for(std::size_t s: util::counter(unpacked_MWR_temp.size())){ + + // If this entry has a unique time them store it for later + if(std::find(MWR_times[dev].begin(), MWR_times[dev].end(), MWR_times_temp[s]) == MWR_times[dev].end()){ + unpacked_MWR[dev].push_back(unpacked_MWR_temp[s]); + MWR_times[dev].push_back(MWR_times_temp[s]); + }//check for unique time + }//Iterate through the unpacked events + }//try + catch (WebAPIException &we) { + //Ignore when we can't find the MWR devices + // they don't always report and the timing of them can be annoying + + }//catch + }// Iterate over all the multiwire devices + }// Iterate over all times + + mf::LogDebug("BNBRetriever") << " Number of MWR[0] times : " << MWR_times[0].size() << std::endl; + mf::LogDebug("BNBRetriever") << " Number of MWR[0]s : " << unpacked_MWR[0].size() << std::endl; + mf::LogDebug("BNBRetriever") << " Number of MWR[1] times : " << MWR_times[1].size() << std::endl; + mf::LogDebug("BNBRetriever") << " Number of MWR[1]s : " << unpacked_MWR[1].size() << std::endl; + mf::LogDebug("BNBRetriever") << " Number of MWR[2] times : " << MWR_times[2].size() << std::endl; + mf::LogDebug("BNBRetriever") << " Number of MWR[2]s : " << unpacked_MWR[2].size() << std::endl; + + return { std::move(MWR_times), std::move(unpacked_MWR) }; +} + + +int sbn::BNBRetriever::matchMultiWireData( + TriggerInfo_t const& triggerInfo, + MWRdata_t const& MWRdata, bool isFirstEventInRun, + std::vector< sbn::BNBSpillInfo >& beamInfos +) const { + + auto const& [ MWR_times, unpacked_MWR ] = MWRdata; // alias + + //Here we will start collecting all the other beamline devices + // First we get the times that the beamline device fired + // we have to pick a specific variable to use + std::vector times_temps = bfp->GetTimeList(fDeviceUsedForTiming); + + // We'll keep track of how many of these spills match to our + // DAQ trigger times + int spill_count = 0; + std::vector matched_MWR; + matched_MWR.resize(3); + + + // Need to handle the first event in a run differently + if(isFirstEventInRun){ + + //We'll remove the spills after our event + int spills_after_our_target = 0; + // iterate through all the spills to find the + // spills that are after our triggered event + for (size_t i = 0; i < times_temps.size(); i++) { + if(times_temps[i] > (triggerInfo.t_current_event+fTimePad)){ + spills_after_our_target++; + } + }//end loop through spill times + + // Remove the spills after our trigger + times_temps.erase(times_temps.end()-spills_after_our_target,times_temps.end()); + + // Remove the spills before the start of our Run + times_temps.erase(times_temps.begin(), times_temps.end() - std::min(int(triggerInfo.number_of_gates_since_previous_event), int(times_temps.size()))); + + }//end fix for "first event" + + // Iterating through each of the beamline times + for (size_t i = 0; i < times_temps.size(); i++) { + + // Only continue if these times are matched to our DAQ time + // plus or minus some time padding, currently using 3.3 ms + // which is half the Booster Rep Rate + + if(!isFirstEventInRun){//We already addressed the "first event" above + if(times_temps[i] > (triggerInfo.t_current_event+fTimePad)){continue;} + if(times_temps[i] <= (triggerInfo.t_previous_event-fTimePad)){continue;} + } + + + //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 = makeBNBSpillInfo(times_temps[i], MWRdata, matched_MWR); + + beamInfos.push_back(std::move(spillInfo)); + + // 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 + + }//end iteration over beam device times + + return spill_count; +} + + +sbn::BNBSpillInfo sbn::BNBRetriever::makeBNBSpillInfo + (double time, MWRdata_t const& MWRdata, std::vector const& matched_MWR) const +{ + + auto const& [ MWR_times, unpacked_MWR ] = MWRdata; // alias + + // initializing all of our device carriers + // device definitions can be found in BNBSpillInfo.h + + double TOR860 = 0; // units e12 protons + double TOR875 = 0; // units e12 protons + double LM875A = 0; // units R/s + double LM875B = 0; // units R/s + double LM875C = 0; // units R/s + double HP875 = 0; // units mm + double VP875 = 0; // units mm + double HPTG1 = 0; // units mm + double VPTG1 = 0; // units mm + double HPTG2 = 0; // units mm + double VPTG2 = 0; // units mm + double BTJT2 = 0; // units Deg C + double THCURR = 0; // units kiloAmps + + double TOR860_time = 0; // units s + + // Here we request all the devices + // since sometimes devices fail to report we'll + // allow each to throw an exception but still move forward + // interpreting these failures will be part of the beam quality analyses + try{bfp->GetNamedData(time, "E:TOR860@",&TOR860,&TOR860_time);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:TOR875",&TOR875);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:LM875A",&LM875A);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:LM875B",&LM875B);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:LM875C",&LM875C);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:HP875",&HP875);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:VP875",&VP875);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:HPTG1",&HPTG1);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:VPTG1",&VPTG1);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:HPTG2",&HPTG2);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:VPTG2",&VPTG2);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:BTJT2",&BTJT2);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + try{bfp->GetNamedData(time, "E:THCURR",&THCURR);}catch (WebAPIException &we) {mf::LogDebug("BNBRetriever")<< "At time : " << time << " " << "got exception: " << we.what() << "\n";} + + //crunch the times + unsigned long int time_closest_int = (int) TOR860_time; + double time_closest_ns = (TOR860_time - time_closest_int)*1e9; + + //Store everything in our data-product + sbn::BNBSpillInfo beamInfo; + beamInfo.TOR860 = TOR860; + beamInfo.TOR875 = TOR875; + beamInfo.LM875A = LM875A; + beamInfo.LM875B = LM875B; + beamInfo.LM875C = LM875C; + beamInfo.HP875 = HP875; + beamInfo.VP875 = VP875; + beamInfo.HPTG1 = HPTG1; + beamInfo.VPTG1 = VPTG1; + beamInfo.HPTG2 = HPTG2; + beamInfo.VPTG2 = VPTG2; + beamInfo.BTJT2 = BTJT2; + beamInfo.THCURR = THCURR; + beamInfo.spill_time_s = time_closest_int; + beamInfo.spill_time_ns = time_closest_ns; + + for(auto const& MWRdata: unpacked_MWR){ + std::ignore = MWRdata; + assert(!MWRdata.empty()); + } + + beamInfo.M875BB = unpacked_MWR[0][matched_MWR[0]]; + beamInfo.M875BB_spill_time_diff = (MWR_times[0][matched_MWR[0]] - time); + beamInfo.M876BB = unpacked_MWR[1][matched_MWR[1]]; + beamInfo.M876BB_spill_time_diff = (MWR_times[1][matched_MWR[1]] - time); + beamInfo.MMBTBB = unpacked_MWR[2][matched_MWR[2]]; + beamInfo.MMBTBB_spill_time_diff = (MWR_times[2][matched_MWR[2]] - time); + + // 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 + + return beamInfo; +} + + +void sbn::BNBRetriever::beginSubRun(art::SubRun& sr) +{ + return; +} + +//____________________________________________________________________________ +void sbn::BNBRetriever::endSubRun(art::SubRun& sr) +{ + // We will add all of the BNBSpillInfo data-products to the + // art::SubRun so it persists + // currently this is ~2.7 kB/event or ~0.07 kB/spill + +mf::LogDebug("BNBRetriever")<< "Total number of DAQ Spills : " << TotalBeamSpills << std::endl; +mf::LogDebug("BNBRetriever")<< "Total number of Selected Spills : " << fOutbeamInfos.size() << std::endl; + + auto p = std::make_unique< std::vector< sbn::BNBSpillInfo > >(fOutbeamInfos); + + sr.put(std::move(p), art::subRunFragment()); + + return; +} + +DEFINE_ART_MODULE(sbn::BNBRetriever) diff --git a/sbncode/BeamSpillInfoRetriever/BNBRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/BNBRetriever/CMakeLists.txt new file mode 100644 index 000000000..99ec10d90 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/BNBRetriever/CMakeLists.txt @@ -0,0 +1,51 @@ +find_ups_product(ifbeam) +find_ups_product(ifbeam_art) + +include_directories( $ENV{IFBEAM_FQ_DIR}/include ) +include_directories( $ENV{LIBWDA_FQ_DIR}/include ) +include_directories( $ENV{IFDHC_FQ_DIR}/inc ) +include_directories( $ENV{IFDH_ART_INC} ) +cet_find_library( IFBEAMSERVICE NAMES IFBeam_service PATHS ENV IFDH_ART_LIB NO_DEFAULT_PATH ) +cet_find_library( IFBEAM NAMES ifbeam PATHS ENV IFBEAM_LIB NO_DEFAULT_PATH ) + +link_libraries( ${LIB_NAME} -L$ENV{BOOST_LIB} -lboost_system ${ROOTLIB} ) + +art_make_library( + LIBRARY_NAME sbn_BNBSpillInfoRetriever_MWRData + SOURCE MWRData.cpp +) + + +simple_plugin(BNBRetriever module + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_ROOT_IO_TFILESERVICE_SERVICE} ${ART_FRAMEWORK_SERVICES} + ${ART_FRAMEWORK_PRINCIPAL} + art_Persistency_Common + art_Utilities canvas + ${MF_MESSAGELOGGER} + ${MF_UTILITIES} + ${FHICLCPP} + cetlib cetlib_except + ${ROOT_EVE_LIB_LIST} + ${ROOT_X3d} + ${ROOT_BASIC_LIB_LIST} + ${Boost_PROGRAM_OPTIONS_LIBRARY} ${Boost_FILESYSTEM_LIBRARY} ${Boost_DATE_TIME_LIBRARY} ${Boost_SERIALIZATION_LIBRARY} + ${Boost_SYSTEM_LIBRARY} + MF_MessageLogger + sbnobj_Common_POTAccounting + #${IFBEAMSERVICE} + #${IFBEAM} + ifbeam + ifdh_art::IFBeam_service + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_Common + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_ICARUS + artdaq_core::artdaq-core_Utilities + sbn_BNBSpillInfoRetriever_MWRData +) + +install_headers() +install_fhicl() +install_source() + diff --git a/sbncode/BeamSpillInfoRetriever/BNBRetriever/MWRData.cpp b/sbncode/BeamSpillInfoRetriever/BNBRetriever/MWRData.cpp new file mode 100644 index 000000000..48d0bda3f --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/BNBRetriever/MWRData.cpp @@ -0,0 +1,71 @@ +#include +#include +#include +#include +#include +#include +#include "MWRData.h" + +using namespace std; + +namespace sbn{ + + std::vector< std::vector < int > > MWRData::unpackMWR(std::string packed_data, std::vector &time_stamp, double timeoffset) const +{ + + std::vector > unpacked_data; + unpacked_data.resize(4); + short data[444]; + + std::vector row(0); + boost::split(row, packed_data, boost::is_any_of(",")); + if (row.size()==447) { + for (int i=3;i<447;i++) { + data[i-3]=atoi(row[i].c_str()); + } + string devname=row[1].substr(0,8); + for (int idev=0;idev<4;idev++) { + mwrpulse_t mwr=getMWRdata(data,idev); + time_stamp.push_back(mwr.sheader.timesec+mwr.sheader.timensec/1000000000.+timeoffset); + for (int ich=0;ich<48;ich++) { + unpacked_data[idev].push_back(mwr.hor[ich]); + } + for (int ich=0;ich<48;ich++) { + unpacked_data[idev].push_back(mwr.ver[ich]); + } + } + } else { + cout <<"BeamSpillInfoRetriever: MRWData: Bad data!"< +namespace sbn{ +class MWRData +{ + typedef struct swicheader_t { + long timesec; + long timensec; + long gpstime1; + long gpstime2; + short boosterevent; + short mievent; + short hz15micnt; + long delta1f; + short pulsemi; + short pulsesc; + } swicheader_t; + + typedef struct mwrpulse_t { + short hor[48]; + short ver[48]; + swicheader_t sheader; + } mwrpulse_t; + + static long flipByte(long data) + { + return ((data>>16)&0x0000FFFF) | ((data<<16)&0xFFFF0000); + } + + mwrpulse_t getMWRdata(short* data, int nblock) const; + + public: + std::vector< std::vector < int > > unpackMWR(std::string packed_data, std::vector &time_stamp, double timeoffset=0) const; +}; +} + +#endif /* #ifndef _MWRDATA_H */ diff --git a/sbncode/BeamSpillInfoRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt new file mode 100644 index 000000000..fba5bc1ff --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/CMakeLists.txt @@ -0,0 +1,7 @@ +add_subdirectory(BNBRetriever) +add_subdirectory(EXTRetriever) +add_subdirectory(job) + +install_headers() +install_fhicl() +install_source() diff --git a/sbncode/BeamSpillInfoRetriever/EXTRetriever/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/EXTRetriever/CMakeLists.txt new file mode 100644 index 000000000..4a3ceaf46 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/EXTRetriever/CMakeLists.txt @@ -0,0 +1,27 @@ + +include_directories( $ENV{IFBEAM_FQ_DIR}/include ) +include_directories( $ENV{LIBWDA_FQ_DIR}/include ) +include_directories( $ENV{IFDHC_FQ_DIR}/inc ) +include_directories( $ENV{IFDH_ART_INC} ) + +simple_plugin(EXTRetriever module + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_FRAMEWORK_SERVICES} + ${ART_FRAMEWORK_PRINCIPAL} + art_Persistency_Common + art_Utilities canvas + ${FHICLCPP} + cetlib cetlib_except + MF_MessageLogger + sbnobj_Common_POTAccounting + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_Common + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays + sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_ICARUS + artdaq_core::artdaq-core_Utilities +) + +install_headers() +install_fhicl() +install_source() + diff --git a/sbncode/BeamSpillInfoRetriever/EXTRetriever/EXTRetriever_module.cc b/sbncode/BeamSpillInfoRetriever/EXTRetriever/EXTRetriever_module.cc new file mode 100644 index 000000000..77d0661e1 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/EXTRetriever/EXTRetriever_module.cc @@ -0,0 +1,155 @@ +//////////////////////////////////////////////////////////////////////// +// Class: EXTRetriever +// Plugin Type: producer +// File: EXTRetriever_module.cc +// +// Created by hand Thurs June 24th 2021 by J. Zennamo (FNAL) +// +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "messagefacility/MessageLogger/MessageLogger.h" +#include "lardata/Utilities/AssociationUtil.h" +#include "lardataobj/Utilities/sparse_vector.h" +#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" +#include "larcorealg/Geometry/Exceptions.h" + +#include "artdaq-core/Data/Fragment.hh" +#include "sbndaq-artdaq-core/Overlays/ICARUS/ICARUSTriggerUDPFragment.hh" + +#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h" +#include "lardata/DetectorInfoServices/DetectorPropertiesService.h" +#include "sbnobj/Common/POTAccounting/EXTCountInfo.h" + +#include +#include +#include +#include +#include +#include + +namespace sbn { + class EXTRetriever; +} + +class sbn::EXTRetriever : public art::EDProducer { +public: + + struct Config { + + using Name = fhicl::Name; + using Comment = fhicl::Comment; + + fhicl::Atom RawDataLabel { + Name{ "raw_data_label" }, + Comment{ "art data product instance name for trigger information (product label is 'daq')" } + }; + + }; // Config + + using Parameters = art::EDProducer::Table; + + + explicit EXTRetriever(Parameters const& params); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + EXTRetriever(EXTRetriever const&) = delete; + EXTRetriever(EXTRetriever&&) = delete; + EXTRetriever& operator=(EXTRetriever const&) = delete; + EXTRetriever& operator=(EXTRetriever&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + void beginSubRun(art::SubRun& sr) override; + void endSubRun(art::SubRun& sr) override; + +private: + std::vector< sbn::EXTCountInfo > fOutExtInfos; + + // input labels + std::string raw_data_label_; + int TotalEXTCounts; + +}; + + +sbn::EXTRetriever::EXTRetriever(Parameters const& params) + : EDProducer{params}, + raw_data_label_(params().RawDataLabel()) +{ + + produces< std::vector< sbn::EXTCountInfo >, art::InSubRun >(); + TotalEXTCounts = 0; +} + +void sbn::EXTRetriever::produce(art::Event& e) +{ + + //Here we read in the artdaq Fragments and extract three pieces of information: + // 1. The time of the current event, t_current_event + // 2. the time of the previously triggered event, t_previous_event (NOTE: Events are non-sequential!) + // 3. the number of beam spills since the previously triggered event, number_of_gates_since_previous_event + + int gate_type = 0; + auto const & raw_data = e.getProduct< std::vector >({ raw_data_label_, "ICARUSTriggerUDP" }); + + unsigned int number_of_gates_since_previous_event = 0; + + for(auto raw_datum : raw_data){ + + icarus::ICARUSTriggerUDPFragment frag(raw_datum); + std::string data = frag.GetDataString(); + char *buffer = const_cast(data.c_str()); + icarus::ICARUSTriggerInfo datastream_info = icarus::parse_ICARUSTriggerString(buffer); + gate_type = datastream_info.gate_type; + number_of_gates_since_previous_event = frag.getDeltaGatesBNB(); + + } + + //We only want to process EXT gates, i.e. type 3 + if(gate_type == 3) + { + // Keep track of the number of beam gates the DAQ thinks + // are in this file + TotalEXTCounts += number_of_gates_since_previous_event; + + //Store everything in our data-product + sbn::EXTCountInfo extInfo; + extInfo.gates_since_last_trigger = number_of_gates_since_previous_event; + + fOutExtInfos.push_back(extInfo); + // 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 + + }//end check on gate type +} //end loop over events + + +void sbn::EXTRetriever::beginSubRun(art::SubRun& sr) +{ + return; +} + +//____________________________________________________________________________ +void sbn::EXTRetriever::endSubRun(art::SubRun& sr) +{ + // We will add all of the EXTCountInfo data-products to the + // art::SubRun so it persists + auto p = std::make_unique< std::vector< sbn::EXTCountInfo > >(fOutExtInfos); + + sr.put(std::move(p)); + + return; +} + +DEFINE_ART_MODULE(sbn::EXTRetriever) diff --git a/sbncode/BeamSpillInfoRetriever/job/CMakeLists.txt b/sbncode/BeamSpillInfoRetriever/job/CMakeLists.txt new file mode 100644 index 000000000..0361ef331 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/CMakeLists.txt @@ -0,0 +1,3 @@ +install_fhicl() +FILE(GLOB fcl_files *.fcl) +install_source( EXTRAS ${fcl_files} ) diff --git a/sbncode/BeamSpillInfoRetriever/job/bnbspillinfo.fcl b/sbncode/BeamSpillInfoRetriever/job/bnbspillinfo.fcl new file mode 100644 index 000000000..4959a0aaa --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/bnbspillinfo.fcl @@ -0,0 +1,16 @@ + +BEGIN_PROLOG + +bnbspillinfo: { + + module_type: "BNBRetriever" + TimePadding: 0.0333 #unit seconds, Booster Rep Rate is 15 Hz, so the closest spill could be 66ms away + URL: "" #keep this blank and we're good + Bundle: "BoosterNeutrinoBeam_read" + MultiWireBundle: "BNBMultiWire" + TimeWindow: "500" #seconds + MWR_TimeWindow: "3601" #seconds + raw_data_label: "daq" + DeviceUsedForTiming: "E:TOR860" +} +END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/job/extcountinfo.fcl b/sbncode/BeamSpillInfoRetriever/job/extcountinfo.fcl new file mode 100644 index 000000000..1ca870996 --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/extcountinfo.fcl @@ -0,0 +1,10 @@ + +BEGIN_PROLOG + +extcountinfo: { + + module_type: "EXTRetriever" + raw_data_label: "daq" + +} +END_PROLOG diff --git a/sbncode/BeamSpillInfoRetriever/job/run_bnbinfo_sbn.fcl b/sbncode/BeamSpillInfoRetriever/job/run_bnbinfo_sbn.fcl new file mode 100644 index 000000000..34f8f689d --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/run_bnbinfo_sbn.fcl @@ -0,0 +1,46 @@ +#include "bnbspillinfo.fcl" + +process_name: BNBInfoGen + +services:{ + + message: { + debugModules: [ "*" ] + destinations: { + LogDebugFile:{ + type: "file" + filename: "debug.log" + append: false + threshold: "DEBUG" + categories: { + default: {} + } + } + } + } + IFBeam:{} +} + + +source: { + +} + +physics: { + producers: { + bnbinfo: @local::bnbspillinfo + } + + simulate: [bnbinfo ] + stream1: [ out1 ] +} + +outputs: { + out1: { + module_type: RootOutput + fileName: "%ifb_%tc_bnbinfo.root" + dataTier: "raw" + compressionLevel: 1 + } +} + diff --git a/sbncode/BeamSpillInfoRetriever/job/run_extinfo_sbn.fcl b/sbncode/BeamSpillInfoRetriever/job/run_extinfo_sbn.fcl new file mode 100644 index 000000000..e55cb29fe --- /dev/null +++ b/sbncode/BeamSpillInfoRetriever/job/run_extinfo_sbn.fcl @@ -0,0 +1,32 @@ +#include "extcountinfo.fcl" + +process_name: EXTInfoGen + +services:{ + IFBeam:{} +} + + +source: { + +} + +physics: { + producers: { + extinfo: @local::extcountinfo + } + + simulate: [extinfo ] + stream1: [ out1 ] + +} + +outputs: { + out1: { + module_type: RootOutput + fileName: "%ifb_%tc_extinfo.root" + dataTier: "raw" + compressionLevel: 1 + } +} + diff --git a/sbncode/CMakeLists.txt b/sbncode/CMakeLists.txt index 64944ac23..0ac05fd7d 100644 --- a/sbncode/CMakeLists.txt +++ b/sbncode/CMakeLists.txt @@ -1,5 +1,7 @@ cet_find_library( BOOST_SERIALIZATION NAMES boost_serialization PATHS ENV BOOST_LIB NO_DEFAULT_PATH ) cet_find_library( BOOST_DATE_TIME NAMES boost_date_time PATHS ENV BOOST_LIB NO_DEFAULT_PATH ) +cet_find_library( IFBEAM NAMES ifbeam PATHS ENV IFBEAM_LIB NO_DEFAULT_PATH ) +cet_find_library( WDA NAMES wda PATHS ENV LIBWDA_LIB NO_DEFAULT_PATH ) find_ups_product(sbnanaobj) @@ -15,6 +17,7 @@ add_subdirectory(Calibration) add_subdirectory(FlashMatch) add_subdirectory(LArRecoProducer) add_subdirectory(TPCReco) +add_subdirectory(BeamSpillInfoRetriever) add_subdirectory(FluxReader) add_subdirectory(EventGenerator) add_subdirectory(PID) diff --git a/ups/product_deps b/ups/product_deps index 9b3da92dc..10559ed6d 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -37,7 +37,7 @@ product version larsoft v09_32_01 sbnobj v09_12_03 sbnanaobj v09_17_01 -sbndaq_artdaq_core v0_07_09of1 +sbndaq_artdaq_core v0_07_09of2 genie_xsec v3_00_04a # list products required ONLY for the build