diff --git a/sbndcode/CRT/CRTReco/CRTStripHitProducer_module.cc b/sbndcode/CRT/CRTReco/CRTStripHitProducer_module.cc index 1f212abd1..1aa7fed0d 100644 --- a/sbndcode/CRT/CRTReco/CRTStripHitProducer_module.cc +++ b/sbndcode/CRT/CRTReco/CRTStripHitProducer_module.cc @@ -23,10 +23,9 @@ #include "sbnobj/SBND/CRT/FEBData.hh" #include "sbnobj/SBND/CRT/CRTStripHit.hh" -#include "sbnobj/SBND/Timing/DAQTimestamp.hh" +#include "sbnobj/SBND/Timing/FrameShiftInfo.hh" #include "sbndcode/Geometry/GeometryWrappers/CRTGeoService.h" -#include "sbndcode/Decoders/PTB/sbndptb.h" #include "sbndcode/Timing/SBNDRawTimingObj.h" #include "sbndcode/ChannelMaps/CRT/CRTChannelMapService.h" @@ -50,10 +49,9 @@ class sbnd::crt::CRTStripHitProducer : public art::EDProducer { void produce(art::Event& e) override; std::vector CreateStripHits(art::Ptr &data, const uint32_t ref_time_s, - const uint32_t ref_time_ns); + const uint32_t ref_time_ns, const uint32_t etrig_time_s, + const uint32_t etrig_time_ns); std::set UnixSet(const std::vector> &datas); - bool SPECTDCReference(art::Event& e, const uint64_t &raw_ts, uint64_t &ref_time); - bool PTBHLTReference(art::Event& e, const uint64_t &raw_ts, uint64_t &ref_time, uint32_t &hlt_code); private: @@ -61,9 +59,13 @@ class sbnd::crt::CRTStripHitProducer : public art::EDProducer { art::ServiceHandle fCRTChannelMapService; std::string fFEBDataModuleLabel; + std::string fFrameShiftModuleLabel; uint16_t fADCThreshold; std::vector fErrorCoeff; bool fAllowFlag1; + bool fApplyETrigWindow; + double fETrigMin; + double fETrigMax; bool fApplyTs0Window; double fTs0Min; double fTs0Max; @@ -72,24 +74,19 @@ class sbnd::crt::CRTStripHitProducer : public art::EDProducer { double fTs1Max; bool fCorrectForDifferentSecond; bool fReferenceTs0; - int fTimingType; - std::string fDAQHeaderModuleLabel; - std::string fDAQHeaderInstanceLabel; - uint32_t fRawTSCorrection; - uint32_t fMaxAllowedRefTimeDiff; - std::string fSPECTDCModuleLabel; - uint32_t fSPECTDCETrigChannel; - std::string fPTBModuleLabel; - std::vector fAllowedPTBHLTs; }; sbnd::crt::CRTStripHitProducer::CRTStripHitProducer(fhicl::ParameterSet const& p) : EDProducer{p} , fFEBDataModuleLabel(p.get("FEBDataModuleLabel")) + , fFrameShiftModuleLabel(p.get("FrameShiftModuleLabel")) , fADCThreshold(p.get("ADCThreshold")) , fErrorCoeff(p.get>("ErrorCoeff")) , fAllowFlag1(p.get("AllowFlag1")) + , fApplyETrigWindow(p.get("ApplyETrigWindow")) + , fETrigMin(p.get("ETrigMin", 0)) + , fETrigMax(p.get("ETrigMax", std::numeric_limits::max())) , fApplyTs0Window(p.get("ApplyTs0Window")) , fTs0Min(p.get("Ts0Min", 0)) , fTs0Max(p.get("Ts0Max", std::numeric_limits::max())) @@ -98,15 +95,6 @@ sbnd::crt::CRTStripHitProducer::CRTStripHitProducer(fhicl::ParameterSet const& p , fTs1Max(p.get("Ts1Max", std::numeric_limits::max())) , fCorrectForDifferentSecond(p.get("CorrectForDifferentSecond")) , fReferenceTs0(p.get("ReferenceTs0")) - , fTimingType(p.get("TimingType", 0)) - , fDAQHeaderModuleLabel(p.get("DAQHeaderModuleLabel", "")) - , fDAQHeaderInstanceLabel(p.get("DAQHeaderInstanceLabel", "")) - , fRawTSCorrection(p.get("RawTSCorrection", 0)) - , fMaxAllowedRefTimeDiff(p.get("MaxAllowedRefTimeDiff", 0)) - , fSPECTDCModuleLabel(p.get("SPECTDCModuleLabel", "")) - , fSPECTDCETrigChannel(p.get("SPECTDCETrigChannel", 4)) - , fPTBModuleLabel(p.get("PTBModuleLabel", "")) - , fAllowedPTBHLTs(p.get>("AllowedPTBHLTs", {})) { produces>(); produces>(); @@ -127,66 +115,39 @@ void sbnd::crt::CRTStripHitProducer::produce(art::Event& e) std::vector> FEBDataVec; art::fill_ptr_vector(FEBDataVec, FEBDataHandle); - uint64_t raw_ts = 0, ref_time = 0; - uint32_t ref_time_s = 0, ref_time_ns = 0; + uint32_t ref_time_s = 0, ref_time_ns = 0, etrig_time_s = 0, etrig_time_ns = 0; if(fReferenceTs0) { - art::Handle DAQHeaderHandle; - e.getByLabel(fDAQHeaderModuleLabel, fDAQHeaderInstanceLabel, DAQHeaderHandle); + art::Handle frameShiftHandle; + e.getByLabel(fFrameShiftModuleLabel, frameShiftHandle); - if(DAQHeaderHandle.isValid()) - { - artdaq::RawEvent rawHeaderEvent = artdaq::RawEvent(*DAQHeaderHandle); - raw_ts = rawHeaderEvent.timestamp() - fRawTSCorrection; - } + if(!frameShiftHandle.isValid()) + throw std::runtime_error("Frame Shift Info object is invalid, check data quality"); - int timingType = fTimingType; - int timingCh = 0; + uint64_t ref_time = frameShiftHandle->FrameDefault(); - if(timingType == 0) - { - uint64_t spec_tdc_ref_time = 0; - - if(SPECTDCReference(e, raw_ts, spec_tdc_ref_time)) - { - ref_time = spec_tdc_ref_time; - timingCh = fSPECTDCETrigChannel; - } - else - ++timingType; - } + ref_time_s = ref_time / sbnd::timing::kSecondInNanoseconds; + ref_time_ns = ref_time % sbnd::timing::kSecondInNanoseconds; - if(timingType == 1) - { - uint64_t ptb_hlt_ref_time = 0; - uint32_t hlt_code = 0; - - if(PTBHLTReference(e, raw_ts, ptb_hlt_ref_time, hlt_code)) - { - ref_time = ptb_hlt_ref_time; - timingCh = hlt_code; - } - else - ++timingType; - } + timingReferenceInfo->timingType = frameShiftHandle->TimingTypeDefault(); + timingReferenceInfo->timingChannel = frameShiftHandle->TimingChannelDefault(); - if(timingType == 2) - { - std::set unix_set = UnixSet(FEBDataVec); - ref_time = unix_set.size() ? *unix_set.rbegin() * static_cast(1e9) : 0; - } + if(timingReferenceInfo->timingType == sbnd::timing::kNoShiftType) + { + std::set unix_set = UnixSet(FEBDataVec); + ref_time_s = unix_set.size() ? *unix_set.rbegin(): 0; + } - ref_time_s = ref_time / static_cast(1e9); - ref_time_ns = ref_time % static_cast(1e9); + uint64_t etrig_time = frameShiftHandle->FrameEtrig(); - timingReferenceInfo->timingType = timingType; - timingReferenceInfo->timingChannel = timingCh; + etrig_time_s = etrig_time / sbnd::timing::kSecondInNanoseconds; + etrig_time_ns = etrig_time % sbnd::timing::kSecondInNanoseconds; } for(auto data : FEBDataVec) { - std::vector newStripHits = CreateStripHits(data, ref_time_s, ref_time_ns); + std::vector newStripHits = CreateStripHits(data, ref_time_s, ref_time_ns, etrig_time_s, etrig_time_ns); for(auto hit : newStripHits) { @@ -203,7 +164,8 @@ void sbnd::crt::CRTStripHitProducer::produce(art::Event& e) } std::vector sbnd::crt::CRTStripHitProducer::CreateStripHits(art::Ptr &data, const uint32_t ref_time_s, - const uint32_t ref_time_ns) + const uint32_t ref_time_ns, const uint32_t etrig_time_s, + const uint32_t etrig_time_ns) { std::vector stripHits; @@ -249,6 +211,27 @@ std::vector sbnd::crt::CRTStripHitProducer::CreateStripH if(fReferenceTs0) t0 -= ref_time_ns; + if(fApplyETrigWindow) + { + double t0_etrig = (int)data->Ts0() + module.t0DelayCorrection; + const int64_t unix_diff = static_cast(ref_time_s) - static_cast(etrig_time_s); + + if(unix_diff < -1 || unix_diff > 1) + { + throw std::runtime_error(Form("Unix timestamps differ by more than 1 (%li)", unix_diff)); + } + + t0_etrig -= etrig_time_ns; + + if(unix_diff == 1) + t0_etrig -= 1e9; + else if(unix_diff == -1) + t0_etrig += 1e9; + + if(t0_etrig < fETrigMin || t0_etrig > fETrigMax) + return stripHits; + } + if(fApplyTs0Window && (t0 < fTs0Min || t0 > fTs0Max)) return stripHits; @@ -309,92 +292,4 @@ std::set sbnd::crt::CRTStripHitProducer::UnixSet(const std::vector> TDCHandle; - e.getByLabel(fSPECTDCModuleLabel, TDCHandle); - - if(!TDCHandle.isValid() || TDCHandle->size() == 0) - return found; - - std::vector> TDCVec; - art::fill_ptr_vector(TDCVec, TDCHandle); - - uint64_t min_diff = std::numeric_limits::max(); - uint64_t min_diff_ts = 0; - - for(auto ts : TDCVec) - { - if(ts->Channel() == fSPECTDCETrigChannel) - { - uint64_t diff = raw_ts > ts->Timestamp() ? raw_ts - ts->Timestamp() : ts->Timestamp() - raw_ts; - - if(diff < min_diff) - { - min_diff = diff; - min_diff_ts = ts->Timestamp(); - found = true; - } - } - } - - if(min_diff > fMaxAllowedRefTimeDiff) - return false; - - ref_time = min_diff_ts; - - return found; -} - -bool sbnd::crt::CRTStripHitProducer::PTBHLTReference(art::Event& e, const uint64_t &raw_ts, uint64_t &ref_time, uint32_t &hlt_code) -{ - bool found = false; - - art::Handle> PTBHandle; - e.getByLabel(fPTBModuleLabel, PTBHandle); - - if(!PTBHandle.isValid() || PTBHandle->size() == 0) - return found; - - std::vector> PTBVec; - art::fill_ptr_vector(PTBVec, PTBHandle); - - uint64_t min_diff = std::numeric_limits::max(); - uint64_t min_diff_ts = 0; - - for(auto ptb : PTBVec) - { - for(auto hlt : ptb->GetHLTriggers()) - { - uint64_t hlt_timestamp = (hlt.timestamp * 20); - std::bitset<32> hlt_word_bitset = std::bitset<32>(hlt.trigger_word); - - for(uint32_t allowed_hlt : fAllowedPTBHLTs) - { - if(hlt_word_bitset[allowed_hlt]) - { - uint64_t diff = raw_ts > hlt_timestamp ? raw_ts - hlt_timestamp : hlt_timestamp - raw_ts; - - if(diff < min_diff) - { - min_diff = diff; - min_diff_ts = hlt_timestamp; - hlt_code = allowed_hlt; - found = true; - } - } - } - } - } - - if(min_diff > fMaxAllowedRefTimeDiff) - return false; - - ref_time = min_diff_ts; - - return found; -} - DEFINE_ART_MODULE(sbnd::crt::CRTStripHitProducer) diff --git a/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl b/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl index 9aa8e2a26..0093419b9 100644 --- a/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl +++ b/sbndcode/CRT/CRTReco/crtrecoproducers_sbnd.fcl @@ -5,6 +5,7 @@ BEGIN_PROLOG crtstriphitproducer_sbnd: { FEBDataModuleLabel: "crtsim" + FrameShiftModuleLabel: "" ADCThreshold: 60 ErrorCoeff: [ 0.26, -0.27, 0.025 ] AllowFlag1: false @@ -17,20 +18,13 @@ crtstriphitproducer_sbnd: crtstriphitproducer_data_sbnd: @local::crtstriphitproducer_sbnd crtstriphitproducer_data_sbnd.FEBDataModuleLabel: "crtdecoder" +crtstriphitproducer_data_sbnd.FrameShiftModuleLabel: "frameshift" crtstriphitproducer_data_sbnd.AllowFlag1: true crtstriphitproducer_data_sbnd.CorrectForDifferentSecond: true crtstriphitproducer_data_sbnd.ReferenceTs0: true -crtstriphitproducer_data_sbnd.TimingType: 0 -crtstriphitproducer_data_sbnd.DAQHeaderModuleLabel: "daq" -crtstriphitproducer_data_sbnd.DAQHeaderInstanceLabel: "RawEventHeader" -crtstriphitproducer_data_sbnd.RawTSCorrection: 367000 -crtstriphitproducer_data_sbnd.MaxAllowedRefTimeDiff: 3000000 -crtstriphitproducer_data_sbnd.SPECTDCModuleLabel: "tdcdecoder" -crtstriphitproducer_data_sbnd.PTBModuleLabel: "ptbdecoder" -crtstriphitproducer_data_sbnd.AllowedPTBHLTs: [ 1, 2, 3, 4, 5, 14, 15 ] -crtstriphitproducer_data_sbnd.ApplyTs0Window: true -crtstriphitproducer_data_sbnd.Ts0Min: -1.71e6 -crtstriphitproducer_data_sbnd.Ts0Max: 1.51e6 +crtstriphitproducer_data_sbnd.ApplyETrigWindow: true +crtstriphitproducer_data_sbnd.ETrigMin: -1.71e6 +crtstriphitproducer_data_sbnd.ETrigMax: 1.51e6 crtclusterproducer_sbnd: { diff --git a/sbndcode/Decoders/CRT/run_crt_ptb_decoder_sbnd.fcl b/sbndcode/Decoders/CRT/run_crt_ptb_decoder_sbnd.fcl deleted file mode 100644 index 7cb1b72b6..000000000 --- a/sbndcode/Decoders/CRT/run_crt_ptb_decoder_sbnd.fcl +++ /dev/null @@ -1,42 +0,0 @@ -#include "crt_decoder_sbnd.fcl" -#include "services_sbnd.fcl" -#include "rootoutput_sbnd.fcl" -#include "crt_channel_map_service_sbnd.fcl" -#include "SBNDPTBDecoderDefaults.fcl" - -process_name: Decoders - -services: -{ - @table::sbnd_services - CRTChannelMapService: @local::crt_channel_map_sbnd -} - -source: -{ - module_type: RootInput -} - -outputs: -{ - out1: - { - @table::sbnd_rootoutput - dataTier: "reconstructed" - } -} - -physics: -{ - producers: - { - crtdecoder: @local::crt_decoder - ptbdecoder: @local::SBNDPTBDecoderDefaults - } - - decoders: [ crtdecoder, ptbdecoder ] - trigger_paths: [ decoders ] - - stream1: [ out1 ] - end_paths: [ stream1 ] -} diff --git a/sbndcode/Decoders/CRT/run_crt_ptb_tdc_decoder_sbnd.fcl b/sbndcode/Decoders/CRT/run_crt_ptb_tdc_decoder_sbnd.fcl index ca58e58b8..8c4f82ba9 100644 --- a/sbndcode/Decoders/CRT/run_crt_ptb_tdc_decoder_sbnd.fcl +++ b/sbndcode/Decoders/CRT/run_crt_ptb_tdc_decoder_sbnd.fcl @@ -4,6 +4,7 @@ #include "crt_channel_map_service_sbnd.fcl" #include "SBNDPTBDecoderDefaults.fcl" #include "spectdc_decoder_sbnd.fcl" +#include "frameshift_sbnd_data.fcl" process_name: Decoders @@ -34,9 +35,10 @@ physics: crtdecoder: @local::crt_decoder ptbdecoder: @local::SBNDPTBDecoderDefaults tdcdecoder: @local::spec_tdc_decoder_sbnd + frameshift: @local::frameshift_data } - decoders: [ crtdecoder, ptbdecoder, tdcdecoder ] + decoders: [ crtdecoder, ptbdecoder, tdcdecoder, frameshift ] trigger_paths: [ decoders ] stream1: [ out1 ] diff --git a/sbndcode/Decoders/CRT/run_crt_ptb_tdc_pmt_decoder_sbnd.fcl b/sbndcode/Decoders/CRT/run_crt_ptb_tdc_pmt_decoder_sbnd.fcl index 6bd07c8b2..1a9e747df 100644 --- a/sbndcode/Decoders/CRT/run_crt_ptb_tdc_pmt_decoder_sbnd.fcl +++ b/sbndcode/Decoders/CRT/run_crt_ptb_tdc_pmt_decoder_sbnd.fcl @@ -5,6 +5,7 @@ #include "SBNDPTBDecoderDefaults.fcl" #include "spectdc_decoder_sbnd.fcl" #include "pmtdecoder.fcl" +#include "frameshift_sbnd_data.fcl" process_name: Decoders @@ -37,9 +38,10 @@ physics: ptbdecoder: @local::SBNDPTBDecoderDefaults tdcdecoder: @local::spec_tdc_decoder_sbnd pmtdecoder: @local::pmtdecoder + frameshift: @local::frameshift_data } - decoders: [ crtdecoder, ptbdecoder, tdcdecoder, pmtdecoder ] + decoders: [ crtdecoder, ptbdecoder, tdcdecoder, frameshift, pmtdecoder ] trigger_paths: [ decoders ] stream1: [ out1 ] diff --git a/sbndcode/Decoders/PMT/CMakeLists.txt b/sbndcode/Decoders/PMT/CMakeLists.txt index 02ea063b9..0b5ec3709 100644 --- a/sbndcode/Decoders/PMT/CMakeLists.txt +++ b/sbndcode/Decoders/PMT/CMakeLists.txt @@ -3,7 +3,7 @@ cet_build_plugin( SBNDPMTDecoder art::module LIBRARIES sbndaq_artdaq_core::sbndaq-artdaq-core_Overlays_SBND sbnobj::SBND_Timing - sbndcode_OpDetSim + sbndcode_OpDetSim lardata::Utilities canvas::canvas lardataobj::RecoBase diff --git a/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc b/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc index 18a56d8d5..62bbc59cd 100644 --- a/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc +++ b/sbndcode/Decoders/PMT/SBNDPMTDecoder_module.cc @@ -29,9 +29,13 @@ #include "sbndaq-artdaq-core/Overlays/Common/CAENV1730Fragment.hh" #include "sbndaq-artdaq-core/Overlays/SBND/PTBFragment.hh" #include "sbndaq-artdaq-core/Overlays/FragmentType.hh" -#include "sbnobj/SBND/Timing/DAQTimestamp.hh" +#include "sbnobj/SBND/Timing/FrameShiftInfo.hh" #include "sbndcode/Decoders/PTB/sbndptb.h" #include "sbndcode/Timing/SBNDRawTimingObj.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" +#include "sbndcode/Calibration/PDSDatabaseInterface/IPMTCalibrationDatabaseService.h" +#include "sbndcode/OpDetSim/TriggerEmulationService.h" +#include "sbndcode/OpDetSim/sbndPDMapAlg.hh" #include "art_root_io/TFileService.h" #include "TH1D.h" @@ -48,8 +52,7 @@ #include #include #include - -#include "sbndcode/OpDetSim/TriggerEmulationService.h" +#include namespace sbndaq { class SBNDPMTDecoder; @@ -78,26 +81,17 @@ class sbndaq::SBNDPMTDecoder : public art::EDProducer { uint32_t get_boardid(artdaq::Fragment & frag); void get_timing(artdaq::Fragment & frag, uint16_t & postpercent, uint32_t & ttt, uint32_t & len, int & tick); + std::vector fill_chmap(sbndDB::PMTCalibrationDatabase const* pmt_calib_db); + private: uint fdebug; - uint ftiming_type; std::vector fcaen_fragment_name; std::string fcaen_module_label; - + std::string fframeshift_module_label; + std::vector fignore_fragid; uint32_t fnominal_length; - - uint32_t fraw_ts_correction; - - std::string fspectdc_product_name; - uint32_t fspectdc_ftrig_ch; - uint32_t fspectdc_etrig_ch; - - std::string fptb_product_name; - std::vector fptb_allowed_hlts; - uint fptb_raw_diff_max; - uint fallowed_time_diff; std::string fpmt_instance_name; @@ -114,6 +108,8 @@ class sbndaq::SBNDPMTDecoder : public art::EDProducer { int fn_maxflashes; uint fn_caenboards; + uint fn_caenpmtchannels; + uint ftiming_caen_offset; uint16_t fthreshold_ftrig; uint16_t fdefault_postpercent; // should be a number between 0 and 100 @@ -122,15 +118,15 @@ class sbndaq::SBNDPMTDecoder : public art::EDProducer { std::vector fset_fragid_map; bool fuse_set_map; - - std::vector fch_map; - int fmon_threshold; // histogram info std::stringstream histname; //raw waveform hist name art::ServiceHandle tfs; uint evt_counter = 0; + + sbndDB::PMTCalibrationDatabase const* fpmt_calib_db; + opdet::sbndPDMapAlg opdetmap; //map for photon detector types }; @@ -141,23 +137,13 @@ sbndaq::SBNDPMTDecoder::SBNDPMTDecoder(fhicl::ParameterSet const& p) fcaen_fragment_name = p.get>("caen_fragment_name"); fcaen_module_label = p.get("caen_module_label","daq"); + fframeshift_module_label = p.get("frameshift_module_label"); fignore_fragid = p.get>("ignore_fragid",{}); fnominal_length = p.get("nominal_length",5000); - - ftiming_type = p.get("timing_type",0); - - fraw_ts_correction = p.get("raw_ts_correction",367000); // ns - - fspectdc_product_name = p.get("spectdc_product_name","tdcdecoder"); - fspectdc_ftrig_ch = p.get("spectdc_ftrig_ch",3); - fspectdc_etrig_ch = p.get("spectdc_etrig_ch",4); - - fptb_product_name = p.get("ptb_product_name","ptbdecoder"); - fptb_allowed_hlts = p.get>("ptb_allowed_hlts",{}); - fptb_raw_diff_max = p.get("ptb_raw_diff_max",3000000); // ns - fallowed_time_diff = p.get("allowed_time_diff",3000); // us!!! + + // configure output labels fpmt_instance_name = p.get("pmtInstanceName","PMTChannels"); fflt_instance_name = p.get("ftrigInstanceName","FTrigChannels"); ftim_instance_name = p.get("timingInstanceName","TimingChannels"); @@ -171,6 +157,8 @@ sbndaq::SBNDPMTDecoder::SBNDPMTDecoder(fhicl::ParameterSet const& p) fn_maxflashes = p.get("n_maxflashes",30); fn_caenboards = p.get("n_caenboards",8); + fn_caenpmtchannels = p.get("n_caenpmtchannels",15); + ftiming_caen_offset = p.get("timing_caen_offset",900); fthreshold_ftrig = p.get("threshold_ftrig",16350); fdefault_postpercent = p.get("default_postpercent",80); ffragid_offset = p.get("fragid_offset",40960); @@ -179,8 +167,7 @@ sbndaq::SBNDPMTDecoder::SBNDPMTDecoder(fhicl::ParameterSet const& p) fset_fragid_map = p.get>("set_fragid_map",{}); fuse_set_map = p.get("use_set_map",false); - fch_map = p.get>("ch_map",{}); - + fpmt_calib_db = lar::providerFrom(); fmon_threshold = p.get("mon_threshold", 15); produces< std::vector< raw::OpDetWaveform > >(fpmt_instance_name); @@ -223,6 +210,8 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) evt_counter++; std::vector> board_frag_v(fn_caenboards); + std::vector ch_map = fill_chmap(fpmt_calib_db); + uint ncont = 0; // counter for number of containers bool found_caen = false; @@ -285,144 +274,22 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) return; } - // create a timing type per event so the default doesn't get overwritten uint64_t event_trigger_time = 0; // in ns - auto timing_type = ftiming_type; - int timing_ch = -1; - - // get the raw event header timestamp - art::Handle header_handle; - uint64_t raw_timestamp = 0; - evt.getByLabel("daq", "RawEventHeader", header_handle); - if ((header_handle.isValid())){ - auto rawheader = artdaq::RawEvent(*header_handle); - raw_timestamp = rawheader.timestamp() - fraw_ts_correction; // includes sec + ns portion - if (fdebug>1) - std::cout << "Raw timestamp (w/ correction) -> " << "ts (ns): " << raw_timestamp % uint64_t(1e9) << ", sec (s): " << raw_timestamp / uint64_t(1e9) << std::endl; - } + + art::Handle frameShiftHandle; + evt.getByLabel(fframeshift_module_label, frameShiftHandle); + + if(!frameShiftHandle.isValid()) + throw std::runtime_error("Frame Shift Info object is invalid, check data quality"); - // get spec tdc product - if (timing_type==0){ - art::Handle> tdcHandle; - evt.getByLabel(fspectdc_product_name,tdcHandle); - bool found_ett = false; - std::vector tdc_etrig_v; - uint64_t min_raw_tdc_diff = uint64_t(1e12); - - if (!tdcHandle.isValid() || tdcHandle->size() == 0){ - if (fdebug>0) std::cout << "No SPECTDC products found." << std::endl; - timing_type++; - } - else{ - if (fdebug>1) std::cout << "SPECTDC (decoded) products found: " << std::endl; - const std::vector tdc_v(*tdcHandle); - - for (size_t i=0; i1)){ - std::cout << " TDC CH " << ch << " -> " - << "name: " << name - << ", ts (ns): " << ts%uint64_t(1e9) - << ", sec (s): " << ts/uint64_t(1e9) - << ", offset: " << offset - << std::endl; - } - if (ch==fspectdc_etrig_ch){ - found_ett = true; - tdc_etrig_v.push_back(ts); - timing_ch = fspectdc_etrig_ch; - } - } - if (found_ett==false) - timing_type++; - else{ - if (tdc_etrig_v.size()==1){ - event_trigger_time = tdc_etrig_v.front()%uint64_t(1e9); - } - else{ // finding the closest ETRIG to the raw timestamp - for (size_t i=0; i < tdc_etrig_v.size(); i++){ - event_trigger_time = tdc_etrig_v.front()%uint64_t(1e9); - auto tdc_etrig = tdc_etrig_v[i]; - uint64_t diff = (tdc_etrig < (raw_timestamp)) ? (raw_timestamp - tdc_etrig) : (tdc_etrig - raw_timestamp); - if (diff < min_raw_tdc_diff){ - event_trigger_time = tdc_etrig%uint64_t(1e9); - min_raw_tdc_diff = diff; - } - } - } - } - } // end else statement for finding tdc - } - if (timing_type==1){ - // get ptb product - art::Handle> ptbHandle; - evt.getByLabel(fptb_product_name,ptbHandle); - bool found_ptb_hlt = false; - if ((!ptbHandle.isValid() || ptbHandle->size() == 0)){ - if (fdebug>0) std::cout << "No PTB products found." << std::endl; - timing_type++; - } - else{ - const std::vector ptb_v(*ptbHandle); - - if (fdebug>1) std::cout << "PTB (decoded) HLTs: " << std::endl; - for (size_t i=0; i ptb_hlt_bitset = std::bitset<32>(trig.trigger_word); - for (auto allowed_hlt : fptb_allowed_hlts){ - if (ptb_hlt_bitset[allowed_hlt]){ - auto diff = (ptb_timestamp < raw_timestamp) ? (raw_timestamp - ptb_timestamp) : (ptb_timestamp - raw_timestamp); - if (diff > fptb_raw_diff_max) continue; - - found_ptb_hlt = true; - event_trigger_time = ptb_timestamp%(uint(1e9)); - timing_ch = allowed_hlt; - if (fdebug>1){ - std::cout << " PTB HLT " << allowed_hlt << " -> " - << "ts (ns): " << ptb_timestamp%(uint(1e9)) - << std::endl; - } - } - } - } - } // end hlt check - } // end of loop over ptb products - if (found_ptb_hlt==false){ - timing_type++; - } - } // end handle validity check - } // end of timing type 1 - if (timing_type==2){ - // CAEN only configuration - // if neither PTB or SPEC products are found, the timestamp will be equal to - // the start of the waveform (according to CAEN TTT and wvfm length) - event_trigger_time = 0; - timing_ch=16; - } + event_trigger_time = frameShiftHandle->FrameDefault()%uint64_t(1e9); + evtTimingInfo->timingType = frameShiftHandle->TimingTypeDefault(); + evtTimingInfo->timingChannel = frameShiftHandle->TimingChannelDefault(); - if (fdebug>1){ - std::cout << "Using " ; - if (timing_type==0) - std::cout << "SPECTDC for event reference time, "; - else if (timing_type==1) - std::cout << "PTB for event reference time, "; - else if (timing_type==2) - std::cout << "CAEN-only for event reference time, "; - std::cout << "channel/word " << timing_ch << ", event trigger time: " << event_trigger_time << " ns" << std::endl; - } - evtTimingInfo->timingType = timing_type; - evtTimingInfo->timingChannel = timing_ch; + // if frameshift has no shift, the timestamp will be equal to + // the start of the waveform (according to CAEN TTT and wvfm length) + if(evtTimingInfo->timingType == sbnd::timing::kNoShiftType) + event_trigger_time = 0; bool extended_flag = false; @@ -561,7 +428,7 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) wvfmHist->SetBinContent(n + 1, (double)combined_wvfm[n]); } double time_diff = (int(iwvfm_start) - int(event_trigger_time))*1e-3; // us - if ((std::abs(time_diff) > fallowed_time_diff) && (timing_type<2)){ + if ((std::abs(time_diff) > fallowed_time_diff) && (event_trigger_time!=0)){ // second rollover between reference time and waveform start if (std::abs(time_diff + 1e6) < fallowed_time_diff) time_diff += 1e6; // us @@ -571,18 +438,16 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) else if (fdebug>1) std::cout << "WARNING: TIME DIFFERENCE IS GREATER THAN " << fallowed_time_diff << " us. Event timestamp: " << event_trigger_time << ", waveform timestamp: " << iwvfm_start << std::endl; } - uint ch; - if (i == 15) - ch = fragid; - else - ch = fch_map.at(fragid*15 + i); - raw::OpDetWaveform waveform(time_diff, ch, combined_wvfm); if ((i == 15) && (foutput_ftrig_wvfm)){ + // for the timing ch of each board, the opdetwaveform channel is equal to the board id + raw::OpDetWaveform waveform(time_diff, fragid, combined_wvfm); fltwvfmVec->push_back(waveform); art::Ptr wvfmPtr = make_fltwvfm_ptr(fltwvfmVec->size()-1); fltTimingAssns->addSingle(brdTimingInfoPtr, wvfmPtr); } else if ((fragid == 8)){ // fyi: this hardcodes the timing caen board index + // for the timing caen, the opdetwaveform channel is offset (900) + chidx + raw::OpDetWaveform waveform(time_diff, ftiming_caen_offset + i, combined_wvfm); if ((foutput_timing_wvfm) && (std::find(fignore_timing_ch.begin(), fignore_timing_ch.end(), i) == fignore_timing_ch.end())){ timwvfmVec->push_back(waveform); art::Ptr wvfmPtr = make_timwvfm_ptr(timwvfmVec->size()-1); @@ -590,6 +455,8 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) } } else{ + // for normal pmt chs, the opdetwaveform channel is from the ch_map + raw::OpDetWaveform waveform(time_diff, ch_map.at(fragid*15+i), combined_wvfm); pmtwvfmVec->push_back(waveform); art::Ptr wvfmPtr = make_pmtwvfm_ptr(pmtwvfmVec->size()-1); pmtTimingAssns->addSingle(brdTimingInfoPtr, wvfmPtr); @@ -615,14 +482,12 @@ void sbndaq::SBNDPMTDecoder::produce(art::Event& evt) } } } // end board loop + board_frag_v.clear(); // loop through flashes art::ServiceHandle tfs; art::ServiceHandle fTriggerService; int PMTPerBoard = fTriggerService->getPMTPerBoard(); - // int fTotalCAENBoards = fTriggerService->getTotalCAENBoards(); - //std::vector< std::vector > MonPulsesAll; - //MonPulseAll.clear(); std::vector MonPulsesFlat; std::vector pulseSizes; MonPulsesFlat.clear(); @@ -812,4 +677,21 @@ uint32_t sbndaq::SBNDPMTDecoder::get_boardid(artdaq::Fragment & frag){ return boardid; } +std::vector sbndaq::SBNDPMTDecoder::fill_chmap(sbndDB::PMTCalibrationDatabase const* pmt_calib_db){ + std::vector ch_map(fn_caenboards*fn_caenpmtchannels,9999); + if (pmt_calib_db==nullptr){ + throw std::runtime_error("PMT Calibration Database pointer is null."); + } + auto nopdets = opdetmap.size(); + for (size_t idet=0; idetgetCAENDigitizer(idet) - 1; + int channel = pmt_calib_db->getCAENDigitizerChannel(idet); + ch_map[board*15 + channel] = idet; + } + return ch_map; +} + DEFINE_ART_MODULE(sbndaq::SBNDPMTDecoder) diff --git a/sbndcode/Decoders/PMT/pmtdecoder.fcl b/sbndcode/Decoders/PMT/pmtdecoder.fcl index ea2f7cc1b..9707dfc6f 100644 --- a/sbndcode/Decoders/PMT/pmtdecoder.fcl +++ b/sbndcode/Decoders/PMT/pmtdecoder.fcl @@ -5,31 +5,16 @@ pmtdecoder: module_type: "SBNDPMTDecoder" debug: 0 # 0=no output, 1=standard, 2=verbose - # output timing system configuration - ## priority system - ## SPEC TDC ETT [0] -> PTB ETT [1] -> CAEN-only [2] - timing_type: 0 # default is SPEC TDC ETT - - # configure CAEN inputs caen_fragment_name: [ "CAENV1730", "ContainerCAENV1730" ] + frameshift_module_label: "frameshift" ignore_fragid: [] # boards to ignore for testing purposes, fragid==8 is the timing board, this is 0-indexed! - nominal_length: 5000 # expected waveform length - raw_ts_correction: 367000 # ns, correction to the rawheader timestamp - - # configure SPEC TDC inputs - spectdc_product_name: "tdcdecoder" # used if the tdc decoder has been run - spectdc_ftrig_ch: 3 # channel number of the spec tdc that stores the flash trigger - spectdc_etrig_ch: 4 # channel number of the spec tdc that stores the event trigger - - # configure PTB inputs; used in case spec etrig is missing - ptb_product_name: "ptbdecoder" - ptb_allowed_hlts: [1,2,3,4,5,14,15] # beam, offbeam, and crossing muon hlts - ptb_trig_raw_diff_max: 3000000 # the maximum allowed time difference in ns + nominal_length: 5000 # ticks, expected waveform length # configure output options # use time difference between reference + CAEN used to determine rollover and/or weirdness allowed_time_diff: 3000 # us!! not nanoseconds + # configure the output data product instance names pmtInstanceName: "PMTChannels" ftrigInstanceName: "FTrigChannels" @@ -41,7 +26,9 @@ pmtdecoder: # more detailed configurables n_maxflashes: 30 # maximum number of flashes anticipated, this is used for memory allocation n_caenboards: 9 # number of caen boards; this is used for memory allocation - # 1-8 = PMT caens, 9th = timing caen + # 1-8 = PMT caens, 9th = timing caen + n_caenpmtchannels: 15 # number of PMT(!) channels per caen board; for memory allocation of chmap + timing_caen_offset: 900 # the channel offset for the timing caen board threshold_ftrig: 16350 # the threshold for the flash trigger channel fragid_offset: 40960 # the offset to subtract to get single digit fragids; if use_set_map is true, should set this to 0 @@ -54,19 +41,6 @@ pmtdecoder: ## to use this, must set `fragid_offset` to 0!!! use_set_map: false set_fragid_map: [3,7,5,9,11,13,15,17,19] # use if `use_set_map` is true: digitizers 0-8 are mapped to these fragids - - # ordered in sets of 15 by increasing fragID - # the PMT_to_Digitizers.txt file in the same directory is the reference for this - # text file is obtained from docdb 33668-v9 - ch_map: [17,16,71,70,41,40,15,14,69,68,13,12,67,66,65, # digitzer 0 - 39,38,11,10,9,8,63,62,37,36,7,6,61,60,64, # digitizer 1 - 95,94,149,148,119,118,93,92,147,146,91,90,145,144,143, # digitizer 2 - 117,116,89,88,87,86,141,140,115,114,85,84,139,138,142, # digitizer 3 - 173,172,227,226,197,196,171,170,225,224,169,168,223,222,167, # digitizer 4 - 195,194,221,220,165,164,219,218,193,192,163,162,217,216,166, # digitizer 5 - 251,250,305,304,275,274,249,248,303,302,247,246,301,300,245, # digitizer 6 - 273,272,299,298,243,242,297,296,271,270,241,240,295,294,244, # digitizer 7 - 900,901,902,903,904,905,906,907,908,909,910,911,912,913,914] # digitizer 8 (TIMING CAEN) } END_PROLOG diff --git a/sbndcode/Decoders/PMT/run_pmtdecoder.fcl b/sbndcode/Decoders/PMT/run_pmtdecoder.fcl index b0c8f5a32..36eca5384 100644 --- a/sbndcode/Decoders/PMT/run_pmtdecoder.fcl +++ b/sbndcode/Decoders/PMT/run_pmtdecoder.fcl @@ -1,6 +1,8 @@ +#include "services_sbnd.fcl" #include "pmtdecoder.fcl" #include "spectdc_decoder_sbnd.fcl" #include "SBNDPTBDecoderDefaults.fcl" +#include "frameshift_sbnd_data.fcl" process_name: PMTDecoder @@ -14,6 +16,7 @@ services: { Baseline: 14250 #Run 1: 14250, Run 2: 14257 MC: false } + @table::sbnd_services # includes PMTCalibDB, TriggerEmulation } #source is a root file @@ -39,6 +42,7 @@ physics: tdcdecoder: @local::spec_tdc_decoder_sbnd pmtdecoder: @local::pmtdecoder ptbdecoder: @local::SBNDPTBDecoderDefaults + frameshift: @local::frameshift_data } filters:{} @@ -46,7 +50,7 @@ physics: analyzers:{} - reco: [tdcdecoder, ptbdecoder, pmtdecoder] + reco: [tdcdecoder, ptbdecoder, frameshift, pmtdecoder] # define the output stream, there could be more than one if using filters diff --git a/sbndcode/Decoders/SPECTDC/run_ptb_tdc_decoder_frameshift.fcl b/sbndcode/Decoders/SPECTDC/run_ptb_tdc_decoder_frameshift.fcl new file mode 100644 index 000000000..e9b6a2050 --- /dev/null +++ b/sbndcode/Decoders/SPECTDC/run_ptb_tdc_decoder_frameshift.fcl @@ -0,0 +1,63 @@ +#include "spectdc_decoder_sbnd.fcl" +#include "SBNDPTBDecoderDefaults.fcl" + +#include "frameshift_sbnd_data.fcl" + +process_name: Decoders + +services: { + TFileService : {fileName: "decoder_hist.root"} + TriggerEmulationService: { + service_type: TriggerEmulationService + MonWidth: 12 + TotalCAENBoards: 8 + PMTPerBoard: 15 + Baseline: 14250 + MC: false + } +} + +#source is a root file +source: +{ + module_type: RootInput + maxEvents: -1 # Number of events to create +} + +outputs: +{ + out1: + { + module_type: RootOutput + fileName: "decoder-frameshift-art.root" + } +} + +physics: +{ + + producers:{ + tdcdecoder: @local::spec_tdc_decoder_sbnd + ptbdecoder: @local::SBNDPTBDecoderDefaults + frameshift: @local::frameshift_data + } + + filters:{} + + analyzers:{} + + + reco: [tdcdecoder, ptbdecoder, frameshift] + + + # define the output stream, there could be more than one if using filters + stream1: [out1] + + # trigger_paths is a keyword and contains the paths that modify the art::event, + # ie filters and producers + trigger_paths: [reco] + + # end_paths is a keyword and contains the paths that do not modify the art::Event, + # ie analyzers and output streams. these all run simultaneously + end_paths: [stream1] +} diff --git a/sbndcode/Decoders/XARAPUCA/SBNDXARAPUCADecoder_module.cc b/sbndcode/Decoders/XARAPUCA/SBNDXARAPUCADecoder_module.cc index ac94c56e1..844b5f6de 100644 --- a/sbndcode/Decoders/XARAPUCA/SBNDXARAPUCADecoder_module.cc +++ b/sbndcode/Decoders/XARAPUCA/SBNDXARAPUCADecoder_module.cc @@ -10,9 +10,10 @@ * @brief Defines and implements the SBNDXARAPUCADecoder class which inherits from an art::EDProducer * as the decoder for V1740B digitizers, intended for the X-ARAPUCAs. * @details The current version of the SBND X-ARAPUCAs decoder implements the updates shown in - * the SBN Document 38475-v1 in the SBN Document Database. + * the SBN Document 43891-v1 in the SBN Document Database. * @note A Python version of the binary decoding is available for testing purposes. You can find * it [here: V1740 binary decoder](https://github.com/aliciavr/V1740_binary_decoder). + * @version 5.0 */ #include "art/Framework/Core/EDProducer.h" @@ -32,12 +33,11 @@ #include "sbndaq-artdaq-core/Overlays/FragmentType.hh" #include "sbndaq-artdaq-core/Overlays/SBND/CAENV1740Fragment.hh" -#include "sbndcode/Decoders/PTB/sbndptb.h" -#include "sbndaq-artdaq-core/Overlays/SBND/PTBFragment.hh" - #include "sbndcode/Timing/SBNDRawTimingObj.h" #include "sbnobj/SBND/Timing/DAQTimestamp.hh" +#include "sbnobj/SBND/Timing/FrameShiftInfo.hh" + #include "lardataobj/RawData/OpDetWaveform.h" #include "art_root_io/TFileService.h" @@ -74,7 +74,6 @@ class sbndaq::SBNDXARAPUCADecoder : public art::EDProducer { void produce(art::Event& e) override; private: - constexpr static uint64_t NANOSEC_IN_SEC = 1'000'000'000; /**< Number of nanoseconds in one second. */ constexpr static uint64_t MICROSEC_IN_NANOSEC = 1'000; /**< Number of nanoseconds in one microsecond. */ constexpr static double NANOSEC_TO_MICROSEC = 1E-3; /**< Conversion factor from nanoseconds to microseconds. */ @@ -87,12 +86,10 @@ class sbndaq::SBNDXARAPUCADecoder : public art::EDProducer { constexpr static uint32_t NUM_CONSECUTIVE_SAMPLES = 3; /**< Number of consecutive samples for the same channel when decoding waveforms. */ constexpr static uint16_t NANOSEC_PER_TICK = 8; /**< Number of nanoseconds per Trigger Time Tag (TTT) tick. */ - constexpr static uint16_t SPEC_TDC_TIMING = 0; /**< Timing reference frame: ETRIG (Event Trigger) timestamp from the SPEC-TDC. */ - constexpr static uint16_t PTB_TIMING = 1; /**< Timing reference frame: HLT (High Level Trigger) timestamp from the PTB (Penn Trigger Board). */ - constexpr static uint16_t CAEN_ONLY_TIMING = 2; /**< Timing reference frame: CAEN-only. */ + constexpr static uint16_t FRAMESHIFT_TIMING = 0; /**< Timing reference frame: FrameShift timestamp. */ + constexpr static uint16_t CAEN_ONLY_TIMING = 1; /**< Timing reference frame: CAEN-only. */ - constexpr static uint16_t HLT_NOT_FOUND = 999; /**< Random value to indicate no HL trigger found. */ - constexpr static uint16_t HLT_TOO_FAR = 1000; /**< Random value to indicate HL trigger found but too far from the raw timestamp. */ + constexpr static uint16_t TTT_DEFAULT = 0; /**< Default integer value for the nominal TTT. */ unsigned int fevent_counter; /**< Event counter. */ @@ -106,48 +103,53 @@ class sbndaq::SBNDXARAPUCADecoder : public art::EDProducer { std::string fwaveforms_instance_name; /**< Name assigned to the product instance containing the waveforms for each channel and board generated by this art::EDProducer. */ std::string ftiming_ref_instance_name; /**< Name assigned to the product instance containing the timing reference information for each event generated by this art::EDProducer. */ - uint16_t ftiming_priority; /**< Timing priority configured: 0 SPEC-TDC RWM, 1 SPEC-TDC ETT, 2 CAEN-only. */ + bool fframeshift_mode; /**< FrameShift mode: 0 use timing priority, 1 use FrameShift timestamp as reference. */ uint16_t factive_timing_frame; /**< Active timing frame while processing each event. */ unsigned int fns_per_sample; /**< Number of nanoseconds per sample. */ - std::string fraw_event_header_product_name; /**< Label identifying the module that created the artdaq::detail::RawEventHeader data product. */ - uint16_t fraw_trig_max_diff; /**< Maximum difference (in ns) allowed between trigger timestamps and raw timestamp. */ - uint32_t fraw_timestamp_correction; /**< Correction (in ns) to be applied to the raw timestamp. */ - - std::string fspectdc_product_name; /**< Label identifying the module containing the SPEC-TDC decoder products. */ - uint32_t fspectdc_ftrig_ch; /**< Channel assigned by the SPEC-TDC to the flash triggers. */ - uint32_t fspectdc_etrig_ch; /**< Channel assigned by the SPEC-TDC to the event triggers. */ - - std::string fptb_product_name; /**< Label identifying the module containing the PTB decoder products. */ - std::vector fallowed_hl_triggers; /**< List of allowed high-level triggers from the PTB. */ + std::string fframeshift_product_name; /**< Label identifying the module containing the FrameShift module products. */ art::ServiceHandle tfs; /**< ServiceHandle object to store the histograms in the decoder_hist.root output file. */ int fstore_debug_waveforms; /**< Number of waveforms to store in the ServiceHandle object for debugging purposes (0: none, -1: all, n: first n waveforms each event). */ - bool fdebug_tdc_handle; /**< If `true` SPEC-TDC information is printed. */ - bool fdebug_ptb_handle; /**< If `true` PTB information is printed. */ + bool fcombine_ext_frag; /**< If `true` combines extended fragments into a single raw::OpDetWaveform object. */ + int32_t fallowed_jittering; /**< Allowed jittering (in ns) between fragments to be combined. */ + + bool fdebug_frameshift_handle; /**< If `true` FrameShift information is printed. */ bool fdebug_fragments_handle; /**< If `true` V1740B CAEN fragments art::Handle information is printed. */ bool fdebug_timing; /**< If `true` timing data is printed. */ bool fdebug_buffer; /**< If `true` the buffer status is printed. */ bool fdebug_waveforms; /**< If `true` waveforms decoding data is printed. */ + bool fdebug_extended_fragments; /**< If `true` extended fragments information is printed. */ + bool fdebug_jittering; /**< If `true` trigger jittering information is printed. */ bool fverbose; /**< If `true` it increases verbosity of console output for detailed processing steps. */ - // Class methods. - void decode_fragment(uint64_t timestamp, std::vector & fragment_indices, const artdaq::Fragment& fragment, std::vector & prod_wvfms); - - bool get_ptb_hlt_timestamp(art::Event& e, uint64_t corr_raw_timestamp, uint64_t & timestamp, uint16_t & hlt_code); - bool get_spec_tdc_etrig_timestamp(art::Event& e, uint64_t corr_raw_timestamp, uint64_t & timestamp); - - void save_prod_wvfm(size_t board_idx, size_t ch, double ini_wvfm_timestamp, const std::vector > & wvfms, std::vector & prod_wvfms); - void save_debug_wvfm(size_t board_idx, size_t fragment_idx, int ch, double ini_wvfm_timestamp, double end_wvfm_timestamp, const std::vector > & wvfms); + // Main processing method. + void decode_fragment(uint64_t timestamp, uint64_t& first_frag_timestamp, size_t& first_frag_idx, int32_t& first_TTT, int32_t& prev_TTT, std::vector & fragment_indices, const artdaq::Fragment& fragment, std::vector & prod_wvfms, std::vector>& wvfms, bool last_one, size_t& appended_samples); + // Timing. + void shift_time(uint64_t TTT_ticks, int64_t TTT_end_ns, uint64_t frag_timestamp, uint64_t timestamp, uint32_t num_samples_per_wvfm, double& ini_wvfm_timestamp, double& end_wvfm_timestamp, size_t& appended_samples); + std::unique_ptr get_frameshift_timestamp(art::Event& e, uint64_t& timestamp); + + // Waveforms decoding. + void decode_waveforms(const artdaq::Fragment& fragment, std::vector>& wvfms, size_t header_size, uint32_t num_channels, uint32_t num_samples_per_wvfm, uint32_t num_words_per_wvfms, uint32_t num_samples_per_group); uint16_t get_sample(uint64_t buffer, uint32_t msb, uint32_t lsb); uint32_t read_word(const uint32_t* & data_ptr); unsigned int get_channel_id(unsigned int board, unsigned int board_channel); - std::string print_timestamp(uint64_t timestamp); - int64_t signed_difference(uint64_t t1, uint64_t t2); + // Combines decoded waveforms. + void append_waveforms(std::vector>& wvfms, const std::vector>& fragment_wvfms, uint32_t num_channels); + + // Dumps and saves waveforms. + void dump_waveforms(std::vector & prod_wvfms, std::vector>& wvfms, size_t first_frag_idx, size_t board_index, uint32_t num_channels, double ini_wvfm_timestamp, double end_wvfm_timestamp); + void save_prod_wvfm(size_t board_idx, size_t ch, double ini_wvfm_timestamp, const std::vector > & wvfms, std::vector & prod_wvfms); + void save_debug_wvfm(size_t board_idx, size_t fragment_idx, int ch, double ini_wvfm_timestamp, double end_wvfm_timestamp, const std::vector > & wvfms); + + // Auxiliary methods. + int64_t signed_difference_UTC_timestamp(uint64_t t1, uint64_t t2); uint64_t abs_difference(uint64_t t1, uint64_t t2); + std::string print_timestamp(uint64_t timestamp); + }; /** @@ -178,33 +180,26 @@ sbndaq::SBNDXARAPUCADecoder::SBNDXARAPUCADecoder(fhicl::ParameterSet const& p) ftiming_ref_instance_name = p.get ("timing_ref_instance_name", ""); // Gets timing configuration. - ftiming_priority = p.get ("timing_priority", SPEC_TDC_TIMING); factive_timing_frame = CAEN_ONLY_TIMING; fns_per_sample = p.get ("ns_per_sample", 16); - // Gets the raw timestamp configuration. - fraw_event_header_product_name = p.get ("raw_timestamp_product_name", "daq"); - fraw_trig_max_diff = p.get ("raw_trig_max_diff", 3000000); - fraw_timestamp_correction = p.get ("raw_timestamp_correction", 367000); - - // SPEC-TDC access configuration. - fspectdc_product_name = p.get("spectdc_product_name", "tdcdecoder"); - fspectdc_ftrig_ch = p.get("spectdc_ftrig_ch", 3); - fspectdc_etrig_ch = p.get("spectdc_etrig_ch", 4); - - // PTB access configuration. - fptb_product_name = p.get("ptb_product_name", "ptbdecoder"); - fallowed_hl_triggers = p.get>("allowed_hl_triggers", {1, 2, 3, 4, 5, 14, 15}); + // FrameShift access configuration. + fframeshift_product_name = p.get("frameshift_product_name", "frameshift"); // Gets the number of waveforms to store in the debug output file. fstore_debug_waveforms = p.get ("store_debug_waveforms", 0); + // Gets the combination of extended fragments option. + fcombine_ext_frag = p.get ("combine_ext_frag", true); + fallowed_jittering = p.get ("allowed_jittering", 64); + // Gets the debug and verbose options. - fdebug_ptb_handle = p.get ("debug_ptb_handle", false); - fdebug_tdc_handle = p.get ("debug_tdc_handle", false); + fdebug_frameshift_handle = p.get ("debug_frameshift_handle", false); fdebug_fragments_handle = p.get ("debug_fragments_handle", false); fdebug_timing = p.get ("debug_timing", false); fdebug_waveforms = p.get ("debug_waveforms", false); + fdebug_extended_fragments = p.get ("debug_extended_fragments", false); + fdebug_jittering = p.get ("debug_jittering", false); fdebug_buffer = p.get ("debug_buffer", false); fverbose = p.get ("verbose", false); @@ -219,13 +214,12 @@ sbndaq::SBNDXARAPUCADecoder::SBNDXARAPUCADecoder(fhicl::ParameterSet const& p) * @details It is the main function of the art::EDProducer module: it processes each event, searches for CAEN V1740 fragments, decodes * them and creates the output products. * 1. It gets the raw timestamp from the artdaq::RawEventHeader product. - * 2. It tries to get a valid timing reference from the selected time frame (SPEC-TDC ETRIG, PTB HLT or CAEN-only). If the selected + * 2. It tries to get a valid timing reference from the selected time frame (FRAMESHIFT or CAEN-only). If the selected * timing frame is not found, it tries to get the next one in the priority list. * 3. It searches for CAEN V1740 fragments in the event, decodes them and creates the output products: a vector of raw::OpDetWaveform. * 4. It dumps the products in the event. */ -void sbndaq::SBNDXARAPUCADecoder::produce(art::Event& e) -{ +void sbndaq::SBNDXARAPUCADecoder::produce(art::Event& e) { if (fverbose) std::cout << "\n > SBNDXARAPUCADecoder::produce: entering the produce function." << std::endl; // Advances the event counter. @@ -233,74 +227,42 @@ void sbndaq::SBNDXARAPUCADecoder::produce(art::Event& e) // Initializes the output instance products. auto prod_wvfms = std::make_unique > (); - auto prod_event_timing_ref_info = std::make_unique (); - prod_event_timing_ref_info->timingType = CAEN_ONLY_TIMING; // CAEN-only timing frame as default. - prod_event_timing_ref_info->timingChannel = 0; + // CAEN-only timing frame as default. + auto prod_event_timing_ref_info = std::make_unique(CAEN_ONLY_TIMING, 0); // CAEN-only timing frame as default. if (fverbose) std::cout << "\n > SBNDXARAPUCADecoder::produce: products initialized." << std::endl; - // Gets the raw timestamp from the artdaq::RawEventHeader product. - uint64_t corr_raw_timestamp = 0; - art::Handle header_handle; - e.getByLabel(fraw_event_header_product_name, "RawEventHeader", header_handle); - if (header_handle.isValid()) { - auto raw_header = artdaq::RawEvent(*header_handle); - corr_raw_timestamp = raw_header.timestamp() - fraw_timestamp_correction; - - if (fverbose | fdebug_tdc_handle | fdebug_ptb_handle | fdebug_timing) std::cout << "\n > SBNDXARAPUCADecoder::produce: RAW timestamp from the event header: " << print_timestamp(corr_raw_timestamp) << "." << std::endl; - if (fdebug_tdc_handle | fdebug_ptb_handle | fdebug_timing) { - std::cout << "\t RAW timestamp: BEFORE correction: " << print_timestamp(raw_header.timestamp()) << " - AFTER correction: " << print_timestamp(corr_raw_timestamp) << "." << std::endl; - } - } - - // [0] ETRIG timestamp, [1] HLT timestamp, [2] CAEN-only timestamp. - uint64_t timestamp = 0; - - // Gets the SPEC-TDC product (if any). - if (ftiming_priority == SPEC_TDC_TIMING) { - if (get_spec_tdc_etrig_timestamp(e, corr_raw_timestamp, timestamp)) { - prod_event_timing_ref_info->timingType = SPEC_TDC_TIMING; - prod_event_timing_ref_info->timingChannel = fspectdc_etrig_ch; - factive_timing_frame = SPEC_TDC_TIMING; - if (fverbose | fdebug_tdc_handle) std::cout << "\n > SBNDXARAPUCADecoder::produce: SPEC-TDC ETRIG reference, " << print_timestamp(timestamp) << ", found close enough to the raw timestamp. Using SPEC-TDC timing frame as reference." << std::endl; - } else { - ftiming_priority = PTB_TIMING; - if (fverbose | fdebug_tdc_handle) std::cout << "\n > SBNDXARAPUCADecoder::produce: SPEC-TDC ETRIG reference not found. Trying to find PTB reference..." << std::endl; - } - } - - // Gets the PTB product (if any). - if (ftiming_priority == PTB_TIMING) { - uint16_t hlt_code = HLT_NOT_FOUND; - if (get_ptb_hlt_timestamp(e, corr_raw_timestamp, timestamp, hlt_code)) { - prod_event_timing_ref_info->timingType = PTB_TIMING; - prod_event_timing_ref_info->timingChannel = hlt_code; - factive_timing_frame = PTB_TIMING; - if (fverbose | fdebug_ptb_handle) std::cout << "\n > SBNDXARAPUCADecoder::produce: PTB HLT reference, " << print_timestamp(timestamp) << ", found close enough to the raw timestamp. Using PTB timing frame as reference." << std::endl; - } else { - ftiming_priority = CAEN_ONLY_TIMING; - if (fverbose | fdebug_ptb_handle) std::cout << "\n > SBNDXARAPUCADecoder::produce: PTB reference not found. Searching for CAEN-only reference..." << std::endl; + uint64_t timestamp = 0; + + // ============= FRAMESHIFT TIMING ============ // + if (auto fs_timing_ref_info = get_frameshift_timestamp(e, timestamp)) { + *prod_event_timing_ref_info = *fs_timing_ref_info; + factive_timing_frame = FRAMESHIFT_TIMING; + if (fverbose | fdebug_frameshift_handle) { + std::cout << "\n > SBNDXARAPUCADecoder::produce: FRAMESHIFT reference, " << print_timestamp(timestamp) << ", found. Using FRAMESHIFT timing frame as reference." << std::endl; + std::cout << "\t FRAMESHIFT reference details: timing type: " << prod_event_timing_ref_info->timingType << ", timing channel: " << prod_event_timing_ref_info->timingChannel << "." << std::endl; } + } else { + if (fverbose | fdebug_frameshift_handle) std::cout << "\n > SBNDXARAPUCADecoder::produce: FrameShift reference not found. Using CAEN-only timing frame as reference." << std::endl; } + // ============================================== // - if (fverbose | fdebug_tdc_handle | fdebug_ptb_handle | fdebug_timing) { - std::cout << "\n > SBNDXARAPUCADecoder::produce: selected timing frame: "; - - if (factive_timing_frame == SPEC_TDC_TIMING) { - std::cout << "SPEC-TDC."; - } else if (factive_timing_frame == PTB_TIMING) { - std::cout << "PTB."; - } else if (ftiming_priority == CAEN_ONLY_TIMING) { - std::cout << "CAEN-only."; - } - - std::cout << " Timestamp: " << print_timestamp(timestamp) <<", timing type: " << prod_event_timing_ref_info->timingType << ", channel: " << prod_event_timing_ref_info->timingChannel << "." << std::endl; - } // Flag to track if valid CAEN fragments are found. bool found_caen = false; - std::vector fragment_indices(fnum_caen_boards, 0); + std::vector fragment_indices(fnum_caen_boards, -1); + std::vector> wvfms; + if (fdebug_extended_fragments) { + std::cout << " Waveforms size: " << wvfms.size() << std::endl; + } + + size_t appended_samples = 0; + size_t first_frag_idx = 0; + int32_t first_TTT = TTT_DEFAULT; + int32_t prev_TTT = 0; + uint64_t first_frag_timestamp = TTT_DEFAULT; + bool last_one = false; if (fverbose | fdebug_fragments_handle) std::cout << "\n > SBNDXARAPUCADecoder::produce: searching for V1740 fragments..." << std::endl; @@ -338,7 +300,8 @@ void sbndaq::SBNDXARAPUCADecoder::produce(art::Event& e) for (size_t f = 0; f < num_caen_fragments; f++) { const artdaq::Fragment fragment = *container_fragment[f].get(); - decode_fragment(timestamp, fragment_indices, fragment, *prod_wvfms); + last_one = f == (num_caen_fragments - 1); + decode_fragment(timestamp, first_frag_timestamp, first_frag_idx, first_TTT, prev_TTT, fragment_indices, fragment, *prod_wvfms, wvfms, last_one, appended_samples); } // End CAEN V1740 fragments loop. } } // End Container fragments loop. @@ -350,7 +313,8 @@ void sbndaq::SBNDXARAPUCADecoder::produce(art::Event& e) // It searches for all CAEN V1740 fragments. for (size_t f = 0; f < frag_handle_size; f++) { const artdaq::Fragment fragment = fragment_handle->at(f); - decode_fragment(timestamp, fragment_indices, fragment, *prod_wvfms); + last_one = f == (frag_handle_size - 1); + decode_fragment(timestamp, first_frag_timestamp, first_frag_idx, first_TTT, prev_TTT, fragment_indices, fragment, *prod_wvfms, wvfms, last_one, appended_samples); } // End CAEN V1740 fragments loop. } } // End extracting CAEN V1740 fragments. @@ -371,232 +335,39 @@ void sbndaq::SBNDXARAPUCADecoder::produce(art::Event& e) } } -/** -* @brief Searches for the SPEC-TDC ETRIG timestamp closest to the raw timestamp if any SPEC-TDC ETRIG product is found in the event. -* @param[in] e The event to be processed. -* @param[in] corr_raw_timestamp The corrected raw timestamp from the artdaq::RawEventHeader product. -* @param[in,out] timestamp The closest ETRIG timestamp to the raw timestamp (if found). -* @return A boolean indicating if a valid ETRIG timestamp was found close enough to the raw timestamp. -* @details It searches for the SPEC-TDC products in the event and looks for the ETRIG timestamps. If any ETRIG -* timestamp is found, it checks which one is the closest to the raw timestamp and if it is close enough (i.e. -* within fraw_trig_max_diff) it returns it as output. -*/ -bool sbndaq::SBNDXARAPUCADecoder::get_spec_tdc_etrig_timestamp(art::Event& e, uint64_t corr_raw_timestamp, uint64_t& timestamp) { - bool ett_found = false; - - art::Handle > tdc_handle; - e.getByLabel(fspectdc_product_name, tdc_handle); - - if (fverbose | fdebug_tdc_handle) std::cout << "\n > SBNDXARAPUCADecoder::get_spec_tdc_etrig_timestamp: searching for SPEC-TDC products..." << std::endl; - - // The art::Handle object is not valid. - if (!tdc_handle.isValid()) { - if (fdebug_tdc_handle) std::cout << "\n\tTDC-SPEC handle not valid for " << fspectdc_product_name << "." << std::endl; - - // The art::Handle object is empty. - } else if (tdc_handle->empty()) { - if (fdebug_tdc_handle) std::cout << "\n\tTDC-SPEC handle is empty." << std::endl; - - // The art::Handle object is valid and not empty. - } else { - if (fdebug_tdc_handle) std::cout << " \n\tDecoded TDC-SPEC products found: " << tdc_handle->size() << " products." << std::endl; - - unsigned int num_event_triggers = 0; - unsigned int num_flash_triggers = 0; - - uint64_t min_diff = std::numeric_limits::max(); - uint64_t closest_tdc_timestamp = 0; - - if (fdebug_tdc_handle) std::cout << "\t\tTDC Channel \t TDC Name \t\t TDC Timestamp [ns] \t\t TDC Offset [ns]"; - - // Loops over the SPEC-TDC products in the event. - for (size_t t = 0; t < tdc_handle->size(); t++) { - const uint64_t tdc_timestamp = tdc_handle->at(t).Timestamp(); // Timestamp of the signal [ns]. - const uint32_t tdc_channel = tdc_handle->at(t).Channel(); // Hardware channel. - - if (fdebug_tdc_handle) { - const std::string tdc_name = tdc_handle->at(t).Name(); // Name of channel input. - const uint64_t tdc_offset = tdc_handle->at(t).Offset(); // Channel specific offset [ns]. - std::cout << "\n\t\t\t " << tdc_channel << " \t\t\t" << tdc_name.c_str() - << "\t\t" << print_timestamp(tdc_timestamp) - << "\t\t " << tdc_offset; - } - - // Gets the closest ETRIG timestamp to the raw timestamp and counts the number of ETRIG and FTRIG. - if (tdc_channel == fspectdc_etrig_ch) { - num_event_triggers++; - - //uint64_t diff = (tdc_timestamp < corr_raw_timestamp) ? (corr_raw_timestamp - tdc_timestamp) : (tdc_timestamp - corr_raw_timestamp); - uint64_t diff = abs_difference(tdc_timestamp, corr_raw_timestamp); - if (diff < min_diff) { - min_diff = diff; - closest_tdc_timestamp = tdc_timestamp; - } - - if (fdebug_tdc_handle) std::cout << "\t\t" << diff << " ns away from the raw timestamp."; - - } else if (tdc_channel == fspectdc_ftrig_ch) { - num_flash_triggers++; - } - - } // End SPEC-TDC products extraction loop. - - if (fdebug_tdc_handle) { - std::cout << "\n\tNumber of event triggers (ETRIG): " << num_event_triggers << "." << std::endl; - std::cout << "\tNumber of flash triggers (FTRIG): " << num_flash_triggers << "." << std::endl; - if (closest_tdc_timestamp) { - std::cout << "\n\tClosest ETRIG timestamp to the raw timestamp: " << print_timestamp(closest_tdc_timestamp) << " with a difference of " << min_diff << " ns." << std::endl; - } else { - std::cout << "\n\tNo ETRIG timestamps found." << std::endl; - } - } - - // Checks if the closest ETRIG timestamp is close enough to the raw timestamp. - if (min_diff > fraw_trig_max_diff) { - closest_tdc_timestamp = 0; - } else { - ett_found = true; - timestamp = closest_tdc_timestamp; - } - } - - tdc_handle.removeProduct(); - return ett_found; -} +// =============== Main fragment processing method =============== // /** - * @brief Searches for the PTB HLT timestamp closest to the raw timestamp if any PTB HLT product is found in the event. - * @param[in] e The event to be processed. - * @param[in] corr_raw_timestamp The corrected raw timestamp from the artdaq::RawEventHeader product. - * @param[in,out] timestamp The closest HLT timestamp to the raw timestamp (if found). - * @param[in,out] hlt_code The code of the closest HL trigger to the raw timestamp (if found). - * @return A boolean indicating if a valid HLT timestamp was found close enough to the raw timestamp. - * @details It searches for the PTB products in the event and looks for the HLT timestamps allowed in the fallowed_hl_triggers - * list. If any HLT timestamp is found, it checks which one is the closest to the raw timestamp and if it is close enough - * (i.e. within fraw_trig_max_diff) it returns it as output along with its trigger code. - */ -bool sbndaq::SBNDXARAPUCADecoder::get_ptb_hlt_timestamp(art::Event& e, uint64_t corr_raw_timestamp, uint64_t& timestamp, uint16_t& hlt_code) { - bool hlt_found = false; - - art::Handle > ptb_handle; - e.getByLabel(fptb_product_name, ptb_handle); - - if (fverbose | fdebug_ptb_handle) std::cout << "\n > SBNDXARAPUCADecoder::get_ptb_hlt_timestamp: PTB timing frame selected. Searching for PTB products..." << std::endl; - - // The art::Handle object is not valid. - if (!ptb_handle.isValid()) { - if (fdebug_ptb_handle) std::cout << "\n\tPTB handle not valid for " << fptb_product_name << "." << std::endl; - - // The art::Handle object is empty. - } else if (ptb_handle->empty()) { - if (fdebug_ptb_handle) std::cout << "\n\tPTB handle is empty." << std::endl; - - // The art::Handle object is valid and not empty. - } else { - if (fdebug_ptb_handle) std::cout << " \n\tDecoded PTB products found: " << ptb_handle->size() << " products."; - - uint64_t min_diff = std::numeric_limits::max(); - uint64_t closest_ptb_timestamp = 0; - uint16_t hlt_type = HLT_NOT_FOUND; - - // Loops over the PTB products in the event. - for (size_t i = 0; i < ptb_handle->size(); i++) { - std::vector hlt_triggers = ptb_handle->at(i).GetHLTriggers(); - - if (!hlt_triggers.empty()) { - if (fdebug_ptb_handle) std::cout << "\n\t\t Number of HL triggers: " << hlt_triggers.size() << "." << std::endl; - - // Loops over the HL triggers in the PTB product. - for (size_t j = 0; j < hlt_triggers.size(); j++) { - raw::ptb::Trigger hlt_trig = hlt_triggers.at(j); - uint64_t raw_ptb_timestamp = hlt_trig.timestamp; - uint64_t ptb_timestamp = raw_ptb_timestamp * 20; - std::bitset<32> hlt_word_bitset = std::bitset<32>(hlt_trig.trigger_word); - - if (fdebug_ptb_handle) { - std::cout << "\t\t\t HLT " << j << ": "; - std::cout << "Timestamp: " << print_timestamp(ptb_timestamp); - std::cout << "\t\tTrigger word: [" << hlt_word_bitset << "]"; - } - - uint32_t allowed_hlt = 0; // Random trigger code. In this case it is used as flag to check if an allowed HLT trigger is found. - // Loops over the allowed HL triggers. - for (size_t k = 0; k < fallowed_hl_triggers.size(); k++) { - if (hlt_word_bitset[fallowed_hl_triggers[k]]) { - allowed_hlt = fallowed_hl_triggers[k]; - uint64_t diff = abs_difference(ptb_timestamp, corr_raw_timestamp); - if (diff < min_diff) { - min_diff = diff; - closest_ptb_timestamp = ptb_timestamp; - hlt_type = allowed_hlt; - } - if (fdebug_ptb_handle) std::cout << " - Allowed HLT: " << allowed_hlt << " - " << diff << " ns away from the raw timestamp." << std::endl; - } - } - if ((fdebug_ptb_handle) && !allowed_hlt) std::cout << std::endl; - } - - } else { - if (fdebug_ptb_handle) std::cout << "\n\t\t No allowed HL triggers found." << std::endl; - } - } // End PTB products extraction loop. - - if (fdebug_ptb_handle) { - if (closest_ptb_timestamp) { - std::cout << "\n\t Closest HLT trigger to the raw timestamp: " << print_timestamp(closest_ptb_timestamp) << " - HLT type: " << hlt_type << ". With a difference of " << min_diff << " ns." << std::endl; - } else { - std::cout << "\n\t No allowed HL triggers found." << std::endl; - } - } - - // Checks if the closest HLT timestamp is close enough to the raw timestamp. - if (min_diff > fraw_trig_max_diff) { - closest_ptb_timestamp = 0; - hlt_code = HLT_TOO_FAR; - } else { - hlt_found = true; - timestamp = closest_ptb_timestamp; - hlt_code = hlt_type; - } - - } - - ptb_handle.removeProduct(); - return hlt_found; -} - -/** - * @brief This function processes a single fragment from a CAEN V1740 and stores the decoded waveforms in the provided product - * container. - * - * @param[in,out] fragment_indices A 1D vector tracking the number of fragments processed for each board. - * @param[in] fragment The input CAEN V1740 fragment containing raw data to be decoded. - * @param[out] prod_wvfms Vector where the decoded waveforms will be stored as raw::OpDetWaveform objects. - * - * @details - * - Identifies the board index corresponding to the fragment ID. - * - Verifies the fragment ID against known boards and ensures it is within a valid range. - * - - * - Decodes the fragment reading raw 32-bit words from the fragment, storing them in a buffer, - * and extracting 12-bit samples. The decoding includes: - * - Accessing metadata for the number of channels and samples. - * - Getting the initial timestamp with respect to the selected timing frame. - * - Binary decoding of the raw waveform. To assign efficiently each sequential sample extracted. These indices formulas are - * applied: - * - The board channel index: - * \f[ - * c = \left( \frac{S}{3} \mod 8 \right) + g \times 8 - * \f] - * - The channel sample index: - * \f[ - * s = (S \mod 3) + \left( \frac{S}{24} \times 3 \right) \mod s_{w}} - * \f] - * Where the group index is computed as \f$ \frac{S}{s_{g}} \f$. - * - Mapping samples to corresponding channels. - * - Populates the output vector (`prod_wvfms`) with decoded waveforms and optionally generates debug waveforms output. + * @brief Decodes a CAEN V1740 fragment, extracts the waveforms and combines them into the output product. + * @param[in] timestamp The valid timestamp used as reference for the event. + * @param[in,out] first_frag_timestamp The nominal timestamp calculated for the fragment being processed. + * @param[in, out] first_frag_idx The first index of the fragments being combined. + * @param[in,out] first_TTT The nominal Trigger Time Tag (TTT) calculated for the fragment being processed. + * @param[in,out] fragment_indices A vector containing the indices of the fragments already processed for each board. + * @param[in] fragment The artdaq::Fragment object to be processed. + * @param[in,out] prod_wvfms The vector of raw::OpDetWaveform objects to be filled with the decoded waveforms. + * @param[in,out] wvfms A 2D vector containing the waveforms for all channels and boards decoded so far. + * @param[in] last_one A boolean flag indicating if the fragment being processed is the last one in the event. + + * @details This method decodes a CAEN V1740 fragment, extracts the waveforms for each channel, shifts them in time according to the + * timing reference and combines them into the output product. + * 1. It checks if the fragment ID corresponds to a valid board. + * 2. It accesses the metadata and header of the fragment to get information about the number of channels, samples, words, etc. + * 3. It calculates the number of samples per waveform and checks if it is a nominal or extended fragment. + * 4. It decodes the waveforms for each channel in the fragment. + * 5. It shifts the waveforms in time according to the timing reference. + * 6. It combines the decoded waveforms into the output product. + * 7. It dumps the waveforms if required. * + * @pre The art::Event object containing the fragment has been processed to get a valid timing reference. + * @post The vector of raw::OpDetWaveform objects is filled with the decoded waveforms from the fragment. + * @see shift_time + * @see decode_waveforms + * @see append_waveforms + * @see dump_waveforms */ -void sbndaq::SBNDXARAPUCADecoder::decode_fragment(uint64_t timestamp, std::vector & fragment_indices, const artdaq::Fragment& fragment, std::vector & prod_wvfms) { + +void sbndaq::SBNDXARAPUCADecoder::decode_fragment(uint64_t timestamp, uint64_t& first_frag_timestamp, size_t& first_frag_idx, int32_t& first_TTT, int32_t& prev_TTT, std::vector & fragment_indices, const artdaq::Fragment& fragment, std::vector & prod_wvfms, std::vector>& wvfms, bool last_one, size_t& appended_samples) { auto fragment_id = fragment.fragmentID() - ffragment_id_offset; auto it = std::find(fboard_id_list.begin(), fboard_id_list.end(), fragment_id); size_t board_idx; @@ -614,6 +385,7 @@ void sbndaq::SBNDXARAPUCADecoder::decode_fragment(uint64_t timestamp, std::vecto } if (valid_fragment) { + fragment_indices[board_idx]++; if (fverbose) std::cout << "\n > SBNDXARAPUCADecoder::decode_fragment: decoding V1740 CAEN fragment " << fragment_indices[board_idx] << " from the board " << board_idx << " (slot " << fboard_id_list[board_idx] << "):" << std::endl; // =============== Accesses Event metadata and Event header for this fragment =============== // @@ -629,7 +401,8 @@ void sbndaq::SBNDXARAPUCADecoder::decode_fragment(uint64_t timestamp, std::vecto // Gets the number of words of the header and the waveforms. uint32_t num_words_per_event = header.eventSize; - uint32_t num_words_per_header = sizeof(CAENV1740EventHeader) / sizeof(uint32_t); + size_t header_size = sizeof(CAENV1740EventHeader); + uint32_t num_words_per_header = header_size / sizeof(uint32_t); uint32_t num_words_per_wvfms = (num_words_per_event - num_words_per_header); uint32_t num_bits_per_all_wvfms = num_words_per_wvfms * BITS_PER_WORD; @@ -637,14 +410,15 @@ void sbndaq::SBNDXARAPUCADecoder::decode_fragment(uint64_t timestamp, std::vecto uint32_t num_remaining_bits = num_bits_per_all_wvfms % BITS_PER_SAMPLE; uint32_t num_samples_per_wvfm = num_samples_per_all_wvfms / num_channels; uint32_t num_samples_per_group = num_samples_per_wvfm * NUM_CHANNELS_PER_GROUP; + uint32_t num_nominal_samples_per_wvfm = metadata->nSamples; if (fverbose | fdebug_waveforms) { - if (metadata->nSamples == num_samples_per_wvfm) { + if (num_nominal_samples_per_wvfm == num_samples_per_wvfm) { std::cout << " > SBNDXARAPUCADecoder::decode_fragment: [NOMINAL FRAGMENT] " << num_samples_per_wvfm << " samples/waveform." << " (" << num_samples_per_group << " samples per group - 8 channels per group -)." << std::endl; } else { std::cout << " > SBNDXARAPUCADecoder::decode_fragment: [EXTENDED FRAGMENT] " << num_samples_per_wvfm << " samples/waveform." << " (" << num_samples_per_group << " samples per group - 8 channels per group -)." << std::endl; } - std::cout << " > SBNDXARAPUCADecoder::decode_fragment: nominal number of samples per waveform: " << metadata->nSamples << "." << std::endl; + std::cout << " > SBNDXARAPUCADecoder::decode_fragment: nominal number of samples per waveform: " << num_nominal_samples_per_wvfm << "." << std::endl; std::cout << " > SBNDXARAPUCADecoder::decode_fragment: number of words for this fragment: " << num_words_per_event << " (Header: " << num_words_per_header << ", Waveform: " << num_words_per_wvfms << ") words." << std::endl; } @@ -658,145 +432,441 @@ void sbndaq::SBNDXARAPUCADecoder::decode_fragment(uint64_t timestamp, std::vecto // =============== Extracts timing information for this fragment =============== // - // Gets the timing information of the CAEN fragment. - int64_t pulse_duration_ns = num_samples_per_wvfm * fns_per_sample; // ns. uint64_t frag_timestamp = fragment.timestamp(); // ns. - int64_t frag_timestamp_s = frag_timestamp / NANOSEC_IN_SEC; // s. - int64_t frag_timestamp_ns = frag_timestamp % NANOSEC_IN_SEC; // ns. - uint32_t TTT_ticks = header.triggerTime(); int64_t TTT_end_ns = TTT_ticks * NANOSEC_PER_TICK; // ns. - // Gets the full TTT timestamp. - uint64_t full_TTT = 0; - // If the fragment timestamp is greater than the TTT end timestamp, it means that rollover occurred. - if (frag_timestamp_ns > TTT_end_ns) { - if (fverbose | fdebug_timing) std::cout << " > SBNDXARAPUCADecoder::decode_fragment: CAEN TTT rollover occurred w.r.t. the fragment timestamp (FTS)." << std::endl; - full_TTT = (frag_timestamp_s + 1) * NANOSEC_IN_SEC + TTT_end_ns; - } else { - full_TTT = frag_timestamp_s * NANOSEC_IN_SEC + TTT_end_ns; + if (fdebug_timing) { + std::cout << "\t\t TTT header.extendedTriggerTime() [TTT_ticks] = " << header.extendedTriggerTime() << " ticks. \t TTT_end_ns = " << print_timestamp(header.extendedTriggerTime() * NANOSEC_PER_TICK) << "." << std::endl; + std::cout << "\t\t TTT header.triggerTimeRollOver(): " << header.triggerTimeRollOver() << std::endl; } - int64_t ref_timestamp = 0; - + //// =============== Start decoding the waveforms =============== // + std::vector > fragment_wvfms(num_channels, std::vector(num_samples_per_wvfm, 0)); + decode_waveforms(fragment, fragment_wvfms, header_size, num_channels, num_samples_per_wvfm, num_words_per_wvfms, num_samples_per_group); + double ini_wvfm_timestamp = 0; double end_wvfm_timestamp = 0; + + if (fverbose) { + std::cout << " > SBNDXARAPUCADecoder::decode_fragment: processing the decoded fragment and combines the extended ones to their nominal ones if needed." << std::endl; + } - // If an ETRIG or HLT timestamp was found it restarts the time from it. Otherwise the CAEN time frame is assigned. - if (factive_timing_frame != CAEN_ONLY_TIMING) { - ref_timestamp = signed_difference(full_TTT, timestamp); // ns. + if (fcombine_ext_frag) { - ini_wvfm_timestamp = (ref_timestamp - pulse_duration_ns) * NANOSEC_TO_MICROSEC; // us. - end_wvfm_timestamp = ref_timestamp * NANOSEC_TO_MICROSEC; // us. - } else { - ref_timestamp = full_TTT; // ns. + // Combination variables and checks. - ini_wvfm_timestamp = ((ref_timestamp - pulse_duration_ns) % NANOSEC_IN_SEC) * NANOSEC_TO_MICROSEC; // ns. - end_wvfm_timestamp = (ref_timestamp % NANOSEC_IN_SEC) * NANOSEC_TO_MICROSEC; // us. - } + bool is_nominal = (num_nominal_samples_per_wvfm == num_samples_per_wvfm); + bool is_first = (fragment_indices[board_idx] == 0); - if (fdebug_timing) { - std::cout << std::fixed << std::setprecision(0); - std::cout << "\t\t ns/tick = " << NANOSEC_PER_TICK << ", ns/sample = " << fns_per_sample << std::endl; - std::cout << "\t\t TTT header.TriggerTime() [TTT_ticks] = " << TTT_ticks << " ticks. \t TTT_end_ns = " << print_timestamp(TTT_end_ns) << "." << std::endl; - std::cout << "\t\t TTT header.extendedTriggerTime() [TTT_ticks] = " << header.extendedTriggerTime() << " ticks. \t TTT_end_ns = " << print_timestamp(header.extendedTriggerTime() * NANOSEC_PER_TICK) << "." << std::endl; - std::cout << "\t\t TTT header.triggerTimeRollOver(): " << header.triggerTimeRollOver() << std::endl; - std::cout << "\t\t Full Fragment timestamp: " << print_timestamp(frag_timestamp) << " = " << frag_timestamp_s << " s " << frag_timestamp_ns << " ns." << std::endl; - std::cout << "\t\t Full TTT - fragment timestamp = "<< abs_difference(full_TTT, frag_timestamp) << " ns." << " Post-percent: " << (double(abs_difference(full_TTT, frag_timestamp)) / double(pulse_duration_ns)) * 100 << "%." << std::endl; - if (factive_timing_frame == SPEC_TDC_TIMING) { - std::cout << "\t ETRIG (SPEC-TDC) timestamp of the fragment: " << std::endl; - std::cout << "\t\t Full UTC ETRIG timestamp: " << print_timestamp(timestamp) << "." << std::endl; - std::cout << "\t\t ETRIG SPEC-TDC difference applied to the CAEN frame (full timestamps): " << print_timestamp(full_TTT) << " - " << print_timestamp(timestamp) << " = " << ref_timestamp << " ns." << std::endl; - } else if (factive_timing_frame == PTB_TIMING) { - std::cout << "\t HLT ETRIG (PTB) timestamp of the fragment: " << std::endl; - std::cout << "\t\t Full UTC HLT ETRIG timestamp: " << print_timestamp(timestamp) << "." << std::endl; - std::cout << "\t\t HLT ETRIG (PTB) difference applied to the CAEN frame (full timestamps): " << print_timestamp(full_TTT) << " - " << print_timestamp(timestamp) << " = " << ref_timestamp << " ns." << std::endl; - } else if (factive_timing_frame == CAEN_ONLY_TIMING) { - std::cout << "\t CAEN trigger timestamp (TTT) of the fragment: " << std::endl; - std::cout << "\t\t Full UTC TTT timestamp: " << print_timestamp(full_TTT) << " = " << full_TTT / NANOSEC_IN_SEC << " s " << TTT_end_ns << " ns." << std::endl; + int32_t frag_length = num_samples_per_wvfm * fns_per_sample; // ns. + + // - Computation of the difference w.r.t. previous TTT difference. + int32_t TTT_dif = 0; + if (prev_TTT > TTT_end_ns) { + if (fverbose | fdebug_timing) std::cout << " > SBNDXARAPUCADecoder::decode_fragment: CAEN TTT rollover occurred." << std::endl; + TTT_dif = prev_TTT - (NANOSEC_IN_SEC - TTT_end_ns); + } else { + TTT_dif = TTT_end_ns - prev_TTT; } - } - if (fverbose | fdebug_timing) { - std::cout << std::fixed << std::setprecision(3) << std::endl; - if (factive_timing_frame == SPEC_TDC_TIMING) { - std::cout << " > SBNDXARAPUCADecoder::decode_fragment: SPEC-TDC time window of " << end_wvfm_timestamp - ini_wvfm_timestamp << " us: [" << ini_wvfm_timestamp << ", " << end_wvfm_timestamp << "] us." << std::endl; - } else if (factive_timing_frame == PTB_TIMING) { - std::cout << " > SBNDXARAPUCADecoder::decode_fragment: PTB time window of " << end_wvfm_timestamp - ini_wvfm_timestamp << " us: [" << ini_wvfm_timestamp << ", " << end_wvfm_timestamp << "] us." << std::endl; - } else { // CAEN_ONLY_TIMING - std::cout << " > SBNDXARAPUCADecoder::decode_fragment: CAEN time window of " << end_wvfm_timestamp - ini_wvfm_timestamp << " us: [" << ini_wvfm_timestamp << ", " << end_wvfm_timestamp << "] us." << std::endl; + // - Computation of the difference between TTTs difference and the length of the fragment in ns (to detect jittering between fragments). + int32_t J = TTT_dif - frag_length; + + bool dump_comb_wvfms = !is_first && (is_nominal || (!is_nominal && (TTT_dif > (frag_length + fallowed_jittering)))); + + if (is_first) { + first_frag_idx = fragment_indices[board_idx]; + first_TTT = TTT_end_ns; + first_frag_timestamp = frag_timestamp; } - } - // =============== Start decoding the waveforms =============== // - if (fverbose) std::cout << " > SBNDXARAPUCADecoder::decode_fragment: binary decoding of the waveforms starting... " << std::endl; - - std::vector > wvfms(num_channels, std::vector(num_samples_per_wvfm, 0)); - - // Absolute sample number [0, TOTAL_NUM_SAMPLES] where TOTAL_NUM_SAMPLES is the total number of samples stored for an event. - uint32_t S = 0; - // Buffer variables. - uint64_t buffer = 0; - uint32_t bits_in_buffer = 0; - - // Data pointer to the beggining of the waveforms stores in the event. - const uint32_t* data_ptr = reinterpret_cast(fragment.dataBeginBytes() + sizeof(CAENV1740EventHeader)); - // Accesses each word, stores it in the buffer and then the samples are extracted from the buffer. - for (size_t j = 0; j < num_words_per_wvfms; j++) { - uint64_t word = read_word(data_ptr); - - // Adds the new word to the buffer and increments the number of bits stored in it. - if (fdebug_buffer) std::cout << buffer << "[word: " << word << "]" << std::endl; - buffer |= word << bits_in_buffer; - bits_in_buffer += BITS_PER_WORD; // bytes * 8 bits/byte - if (fdebug_buffer) std::cout << " +" << buffer << " [bits in buffer: "<< bits_in_buffer << "]" << std::endl; - - // Obtains 12-bit sequences from the buffer and assigns each sample to the channel and channel sample it belongs to. - while (bits_in_buffer >= BITS_PER_SAMPLE) { - // Computes board channel, channel sample and group channel and assigns the sample to those indices. - uint32_t g = (S / num_samples_per_group); // Group index. - uint32_t c = ((S / NUM_CONSECUTIVE_SAMPLES) % NUM_CHANNELS_PER_GROUP) + g * NUM_GROUPS; // Channel index. - uint32_t s = (S % NUM_CONSECUTIVE_SAMPLES) + ((S / NUM_SAMPLES_PER_ROUND) * NUM_CONSECUTIVE_SAMPLES) % num_samples_per_wvfm; // Sample/channel index. - uint16_t sample = get_sample(buffer, BITS_PER_SAMPLE - 1, 0); - wvfms[c][s] = sample; - if (fdebug_waveforms) std::cout << "\tSample: " << sample << "\tg: " << g << "\tch: " << c << "\ts:" << s << "\tS: " << S << std::endl; - - // Updates the buffer status removing the read bits and decreasing the number of bits stored in it. - buffer >>= BITS_PER_SAMPLE; - bits_in_buffer -= BITS_PER_SAMPLE; - if (fdebug_buffer) std::cout << " -" << buffer << " [bits in buffer: "<< bits_in_buffer << "]" << std::endl; - - // Increments the absolute sample step. - S++; + if (dump_comb_wvfms) { + shift_time(TTT_ticks, first_TTT, first_frag_timestamp, timestamp, num_nominal_samples_per_wvfm, ini_wvfm_timestamp, end_wvfm_timestamp, appended_samples); + dump_waveforms(prod_wvfms, wvfms, first_frag_idx, board_idx, num_channels, ini_wvfm_timestamp, end_wvfm_timestamp); + first_frag_idx = fragment_indices[board_idx]; + first_TTT = TTT_end_ns; + first_frag_timestamp = frag_timestamp; + if (fdebug_jittering) std::cout << "NC - "; + } else { + if (fdebug_jittering) std::cout << "C - "; } - } + append_waveforms(wvfms, fragment_wvfms, num_channels); + appended_samples = wvfms[0].size() - num_nominal_samples_per_wvfm; + + if (last_one) { + if (fdebug_extended_fragments) std::cout << "\t\t LAST fragment " << std::endl; + shift_time(TTT_ticks, first_TTT, first_frag_timestamp, timestamp, num_nominal_samples_per_wvfm, ini_wvfm_timestamp, end_wvfm_timestamp, appended_samples); + dump_waveforms(prod_wvfms, wvfms, first_frag_idx, board_idx, num_channels, ini_wvfm_timestamp, end_wvfm_timestamp); + } + + if (fdebug_jittering) { + if (is_first) { + std::cout << "\t\t FIRST "; + } + if (is_nominal) { + int32_t TTT_ini_ns = TTT_end_ns - static_cast(num_nominal_samples_per_wvfm * fns_per_sample); + std::cout << "NOMINAL fragment " << fboard_id_list[board_idx] << " - " << fragment_indices[board_idx] << " (" << num_samples_per_wvfm << "): [" << TTT_ini_ns << ", " << TTT_end_ns << "] diff with prev: " << TTT_dif << " ns, length: " << frag_length << " ns, diff_with_length (jitt): " << J << std::endl; + } else { + int32_t TTT_ini_ns = TTT_end_ns - static_cast(num_nominal_samples_per_wvfm * fns_per_sample); + std::cout << "EXTENDED fragment " << fboard_id_list[board_idx] << " - " << fragment_indices[board_idx] << " (" << num_samples_per_wvfm << "): [" << TTT_ini_ns << ", " << TTT_end_ns << "] diff with prev: " << TTT_dif << " ns, length: " << frag_length << " ns, diff_with_length (jitt): " << J << std::endl; + } + } + + // WARNINGS + if (TTT_dif < frag_length - fallowed_jittering) { + if (!is_first) { + if (!is_nominal) std::cout << "\t\t WARNING: EXTENDED fragment too close! This fragment is being combined." << std::endl; + if (is_nominal) std::cout << "\t\t WARNING: NOMINAL fragment too close! This fragment is not being combined." << std::endl; + } + } + + prev_TTT = TTT_end_ns; - // The decoded waveforms are dumped into two products: - // - A xarapucadecoder-art.root file with the OpDetWaveforms as the product of this producer for further analysis. - // - A decoder_hist.root file gathering a waveform histograms. - if (fverbose) std::cout << " > SBNDXARAPUCADecoder::decode_fragment: binary decoding complete, dumping products..." << std::endl; + } else { // Combination of extended fragments disabled. + first_frag_idx = fragment_indices[board_idx]; + shift_time(TTT_ticks, TTT_end_ns, frag_timestamp, timestamp, num_samples_per_wvfm, ini_wvfm_timestamp, end_wvfm_timestamp, appended_samples); + dump_waveforms(prod_wvfms, fragment_wvfms, first_frag_idx, board_idx, num_channels, ini_wvfm_timestamp, end_wvfm_timestamp); + } + } + +} + +// =============== Timing functions =============== // + +/** + * @brief This function shifts the initial and end timestamps of a waveform based on the provided timing information. + * + * @param[in] TTT_ticks The trigger time tag in ticks from the CAEN V1740 header. + * @param[in] TTT_end_ns The end time of the Trigger Time Tag (TTT) in nanoseconds. + * @param[in] frag_timestamp The timestamp of the fragment in nanoseconds. + * @param[in] timestamp The reference timestamp (from FRAMESHIFT or CAEN-only) in nanoseconds. + * @param[in] num_samples_per_wvfm The number of samples per waveform. + * @param[out] ini_wvfm_timestamp The initial timestamp of the waveform in microseconds (output). + * @param[out] end_wvfm_timestamp The end timestamp of the waveform in microseconds (output). + * + * @details + * This function calculates the initial and end timestamps for a waveform based on the provided timing information. It takes into account + * potential rollovers in the trigger time tag and adjusts the timestamps accordingly. The function supports different timing frames, including + * FRAMESHIFT and CAEN-only timing. The calculated timestamps are returned in microseconds. + * + */ +void sbndaq::SBNDXARAPUCADecoder::shift_time(uint64_t TTT_ticks, int64_t TTT_end_ns, uint64_t frag_timestamp, uint64_t timestamp, uint32_t num_samples_per_wvfm, double& ini_wvfm_timestamp, double& end_wvfm_timestamp, size_t& appended_samples) { + + int64_t pulse_duration_ns = num_samples_per_wvfm * fns_per_sample; // ns. + int64_t frag_timestamp_s = frag_timestamp / NANOSEC_IN_SEC; // s. + int64_t frag_timestamp_ns = frag_timestamp % NANOSEC_IN_SEC; // ns. + + // Gets the UNIX TTT timestamp. + uint64_t TTT_UNIX_timestamp = 0; + // If the fragment timestamp is greater than the TTT end timestamp, it means that rollover occurred. + if (frag_timestamp_ns > TTT_end_ns) { + if (fverbose | fdebug_timing) std::cout << " > SBNDXARAPUCADecoder::decode_fragment: CAEN TTT rollover occurred w.r.t. the fragment timestamp (FTS)." << std::endl; + TTT_UNIX_timestamp = (frag_timestamp_s + 1) * NANOSEC_IN_SEC + TTT_end_ns; + } else { + TTT_UNIX_timestamp = frag_timestamp_s * NANOSEC_IN_SEC + TTT_end_ns; + } + + int64_t ref_timestamp = 0; + int64_t appended_ns = static_cast(appended_samples) * static_cast(fns_per_sample); // this is just for the debug output. + + // If an frameshift timestamp was found it restarts the time from it. Otherwise the CAEN time frame is assigned. + if (factive_timing_frame != CAEN_ONLY_TIMING) { + ref_timestamp = signed_difference_UTC_timestamp(TTT_UNIX_timestamp, timestamp); // ns. + ini_wvfm_timestamp = (ref_timestamp - pulse_duration_ns) * NANOSEC_TO_MICROSEC; // us. + end_wvfm_timestamp = (ref_timestamp + appended_ns) * NANOSEC_TO_MICROSEC; // us. + } else { + ref_timestamp = TTT_UNIX_timestamp; // ns. + + ini_wvfm_timestamp = ((ref_timestamp - pulse_duration_ns) % NANOSEC_IN_SEC) * NANOSEC_TO_MICROSEC; // ns. + end_wvfm_timestamp = ((ref_timestamp + appended_ns)% NANOSEC_IN_SEC) * NANOSEC_TO_MICROSEC; // us. + } + appended_samples = 0; + if (fdebug_timing) { + std::cout << std::fixed << std::setprecision(0); + std::cout << "\t\t ns/tick = " << NANOSEC_PER_TICK << ", ns/sample = " << fns_per_sample << std::endl; + std::cout << "\t\t TTT header.TriggerTime() [TTT_ticks] = " << TTT_ticks << " ticks. \t TTT_end_ns = " << print_timestamp(TTT_end_ns) << "." << std::endl; + std::cout << "\t\t Full Fragment timestamp: " << print_timestamp(frag_timestamp) << " = " << frag_timestamp_s << " s " << frag_timestamp_ns << " ns." << std::endl; + std::cout << "\t\t Full TTT - fragment timestamp = "<< abs_difference(TTT_UNIX_timestamp, frag_timestamp) << " ns." << " Post-percent: " << (double(abs_difference(TTT_UNIX_timestamp, frag_timestamp)) / double(pulse_duration_ns)) * 100 << "%." << std::endl; + if (factive_timing_frame == FRAMESHIFT_TIMING) { + std::cout << "\t FRAMESHIFT timestamp of the fragment: " << std::endl; + std::cout << "\t\t Full UTC FRAMESHIFT timestamp: " << print_timestamp(timestamp) << "." << std::endl; + std::cout << "\t\t FRAMESHIFT difference applied to the CAEN frame (UNIX timestamps): " << print_timestamp(TTT_UNIX_timestamp) << " - " << print_timestamp(timestamp) << " = " << ref_timestamp << " ns." << std::endl; + } else { + std::cout << "\t CAEN trigger timestamp (TTT) of the fragment: " << std::endl; + std::cout << "\t\t Full UTC TTT timestamp: " << print_timestamp(TTT_UNIX_timestamp) << " = " << TTT_UNIX_timestamp / NANOSEC_IN_SEC << " s " << TTT_end_ns << " ns." << std::endl; + } + } + + if (fverbose | fdebug_timing) { + std::cout << std::fixed << std::setprecision(3); + if (factive_timing_frame == FRAMESHIFT_TIMING) { + std::cout << " > SBNDXARAPUCADecoder::decode_fragment: FRAMESHIFT time window of " << end_wvfm_timestamp - ini_wvfm_timestamp << " us: [" << ini_wvfm_timestamp << ", " << end_wvfm_timestamp << "] us." << std::endl; + } else { // CAEN_ONLY_TIMING + std::cout << " > SBNDXARAPUCADecoder::decode_fragment: CAEN time window of " << end_wvfm_timestamp - ini_wvfm_timestamp << " us: [" << ini_wvfm_timestamp << ", " << end_wvfm_timestamp << "] us." << std::endl; + } + std::cout << " > SBNDXARAPUCADecoder::decode_fragment: TTT_end_ticks = " << TTT_ticks << " ticks. \t TTT_end_ns = " << print_timestamp(TTT_end_ns) << "." << std::endl; + } +} + +/** + * @brief Searches for the frameshift timestamp in the event if any frameshift product is found and returns it as a TimingReferenceInfo object. + * @param[in] e The event to be processed. + * @param[in,out] timestamp The frameshift timestamp (if found). + * @param[in,out] timing_type The type of the frameshift timestamp (if found). + * @param[in,out] timing_channel The channel of the frameshift timestamp (if found). + * @return A unique pointer to a TimingReferenceInfo object containing the frameshift timestamp information (if found), or nullptr if no valid frameshift product is found. + * @details It searches for the frameshift products in the event and looks for the default frameshift timestamp, type and channel. If any valid frameshift product is found, it returns its information as a TimingReferenceInfo object. + */ +std::unique_ptr sbndaq::SBNDXARAPUCADecoder::get_frameshift_timestamp(art::Event& e, uint64_t& timestamp) { + std::unique_ptr fs_timing_ref_info = nullptr; + + art::Handle frameshiftHandle; + e.getByLabel(fframeshift_product_name, frameshiftHandle); + + if (fverbose || fdebug_frameshift_handle) std::cout << "\n > SBNDXARAPUCADecoder::get_frameshift_timestamp: Frameshift timing frame selected. Searching for Frameshift products..." << std::endl; + + if (frameshiftHandle.isValid()) { - uint32_t num_debug_wvfms; + uint16_t timing_type = frameshiftHandle->TimingTypeDefault(); + uint16_t timing_channel = frameshiftHandle->TimingChannelDefault(); - if (fstore_debug_waveforms == -1) { - num_debug_wvfms = num_channels; + if (timing_type != sbnd::timing::kNoShiftType) { // Only if timing type is valid, the timestamp is considered valid. + timestamp = frameshiftHandle->FrameDefault(); + fs_timing_ref_info = std::make_unique(timing_type, timing_channel); } else { - num_debug_wvfms = std::min(num_channels, fstore_debug_waveforms); + timestamp = 0; } + + } else { + if (fdebug_frameshift_handle) std::cout << "\n\tFrameshift handle not valid for " << fframeshift_product_name << "." << std::endl; + throw std::runtime_error("XA Decoder: Frame Shift Info object is invalid, check data quality"); + } + + return fs_timing_ref_info; +} - uint32_t ch; +// =============== Decodes the waveforms =============== // - for (ch = 0; ch < num_debug_wvfms; ch++) { - save_prod_wvfm(board_idx, ch, ini_wvfm_timestamp, wvfms, prod_wvfms); - save_debug_wvfm(board_idx, fragment_indices[board_idx], ch, ini_wvfm_timestamp, end_wvfm_timestamp, wvfms); +/** + * @brief Decodes the waveforms from a CAEN V1740 fragment (binary decoding stage). + * + * @param[in] fragment The input CAEN V1740 fragment containing raw data to be decoded. + * @param[out] wvfms A 2D vector where the decoded waveforms will be stored. Each inner vector corresponds to a channel's waveform. + * @param[in] header_size The size of the event header in bytes. + * @param[in] num_channels The number of channels in the fragment. + * @param[in] num_samples_per_wvfm The number of samples per waveform for each channel. + * @param[in] num_words_per_wvfms The total number of 32-bit words containing waveform data in the fragment. + * @param[in] num_samples_per_group The number of samples per group of channels (8 channels per group). + * + * @details This function reads raw 32-bit words from the fragment, storing them in a buffer, and extracts 12-bit samples. + * It assigns each sample to the appropriate channel and sample index based on its position in the sequence. The decoding process includes: + * - Initializing a buffer to hold incoming data and tracking the number of bits currently stored. + * - Iterating over each 32-bit word in the waveform data section of the fragment, adding it to the buffer. + * - Extracting 12-bit samples from the buffer as long as there are enough bits available, and assigning them to their respective channels and sample indices using calculated formulas. + * The function ensures that all samples are correctly mapped to their channels, taking into account the grouping of channels and the interleaving of samples. + * @note The function assumes that the input fragment is valid and contains the expected structure for a CAEN V1740 device. + * @see SBN Document 38475-v1 for more details on the binary decoding. + */ +void sbndaq::SBNDXARAPUCADecoder::decode_waveforms(const artdaq::Fragment& fragment, std::vector>& wvfms, size_t header_size, uint32_t num_channels, uint32_t num_samples_per_wvfm, uint32_t num_words_per_wvfms, uint32_t num_samples_per_group) { + // =============== Start decoding the waveforms =============== // + if (fverbose) std::cout << " > SBNDXARAPUCADecoder::decode_fragment: binary decoding of the waveforms starting... " << std::endl; + + // Absolute sample number [0, TOTAL_NUM_SAMPLES] where TOTAL_NUM_SAMPLES is the total number of samples stored for an event. + uint32_t S = 0; + // Buffer variables. + uint64_t buffer = 0; + uint32_t bits_in_buffer = 0; + + // Data pointer to the beggining of the waveforms stores in the event. + const uint32_t* data_ptr = reinterpret_cast(fragment.dataBeginBytes() + header_size); + // Accesses each word, stores it in the buffer and then the samples are extracted from the buffer. + for (size_t j = 0; j < num_words_per_wvfms; j++) { + uint64_t word = read_word(data_ptr); + + // Adds the new word to the buffer and increments the number of bits stored in it. + if (fdebug_buffer) std::cout << buffer << "[word: " << word << "]" << std::endl; + buffer |= word << bits_in_buffer; + bits_in_buffer += BITS_PER_WORD; // bytes * 8 bits/byte + if (fdebug_buffer) std::cout << " +" << buffer << " [bits in buffer: "<< bits_in_buffer << "]" << std::endl; + + // Obtains 12-bit sequences from the buffer and assigns each sample to the channel and channel sample it belongs to. + while (bits_in_buffer >= BITS_PER_SAMPLE) { + // Computes board channel, channel sample and group channel and assigns the sample to those indices. + uint32_t g = (S / num_samples_per_group); // Group index. + uint32_t c = ((S / NUM_CONSECUTIVE_SAMPLES) % NUM_CHANNELS_PER_GROUP) + g * NUM_GROUPS; // Channel index. + uint32_t s = (S % NUM_CONSECUTIVE_SAMPLES) + ((S / NUM_SAMPLES_PER_ROUND) * NUM_CONSECUTIVE_SAMPLES) % num_samples_per_wvfm; // Sample/channel index. + uint16_t sample = get_sample(buffer, BITS_PER_SAMPLE - 1, 0); + wvfms[c][s] = sample; + if (fdebug_waveforms) std::cout << "\tSample: " << sample << "\tg: " << g << "\tch: " << c << "\ts:" << s << "\tS: " << S << std::endl; + + // Updates the buffer status removing the read bits and decreasing the number of bits stored in it. + buffer >>= BITS_PER_SAMPLE; + bits_in_buffer -= BITS_PER_SAMPLE; + if (fdebug_buffer) std::cout << " -" << buffer << " [bits in buffer: "<< bits_in_buffer << "]" << std::endl; + + // Increments the absolute sample step. + S++; } + } +} + +/** + * @brief Extract a sample from a 64-bit buffer using the specified bit positions. + * + * @param[in] buffer An unsigned 64-bit integer which represents a temporal buffer for the read words and where the samples are extracted from. + * @param[in] msb An unsigned 32-bit integer representing the most significative bit (MSB) where the readout from the buffer paramter. + * @param[in] lsb An unsigned 32-bit integer representing the less significative bit (LSB) from we end read + * + * @details The function shifts the buffer to the right by the number of positions specified by `lsb` so that the least significant bit of the + * sample aligns with bit 0. It then applies a mask to isolate the bits between `lsb` and `msb`, inclusive. + * + * @return The extracted sample as a 16-bit unsigned integer. + */ +uint16_t sbndaq::SBNDXARAPUCADecoder::get_sample(uint64_t buffer, uint32_t msb, uint32_t lsb) { + uint64_t mask = (1U << (msb - lsb + 1)) - 1; + uint64_t sample = buffer >> lsb; + return sample & mask; +} - for (;ch < num_channels; ch++) { - save_prod_wvfm(board_idx, ch, ini_wvfm_timestamp, wvfms, prod_wvfms); +/** + * @brief Read a 32-bit word from the data pointer and advances the pointer. + * + * @param[in, out] data_ptr A reference to a pointer pointing to the current position in the data. + * + * @details This function retrieves a 32-bit word from the memory location pointed to by `data_ptr`. After reading, it advances `data_ptr` to + * the next 32-bit word location. + * + * @return The 32-bit word read from the location pointed to by `data_ptr`. + */ +uint32_t sbndaq::SBNDXARAPUCADecoder::read_word(const uint32_t* & data_ptr) { + uint32_t word = *data_ptr; + data_ptr += 1; + return word; +} + +/** + * @brief Generates a unique global channel identifier using the board slot and the channel number of that board. + * + * @param[in] board Index of the board in `fboard_id_list` from which to derive the board slot. + * @param[in] board_channel The specific channel number on the given board. + * + * @details This function computes a `channel_id` by combining the board slot and the specific + * channel number on that board. + * The unique identifier `channel_id` (\f$ CH_{ID} $\f) is computed as follows: + * \f[ + * CH\_{ID} = B\_{ID} \times 100 + CH\_{B} + * \f] + * + * Where: + * - \f$ B\_{ID} \f$ is the fragment ID retrieved from `fboard_id_list` based on the slot. + * - \f$ CH\_B \f$ is the channel number on that board. + * + * @return A unique identifier for the specified channel as an unsigned integer. + */ +unsigned int sbndaq::SBNDXARAPUCADecoder::get_channel_id(unsigned int board, unsigned int board_channel) { + unsigned int channel_id = fboard_id_list[board] * 100 + board_channel; + return channel_id; +} + +// =============== Combines the waveforms from extended fragments =============== // + +/** + * @brief Combines the waveforms from the current fragment with the previously stored waveforms. + * @param[in,out] wvfms A 2D vector containing the (combined if it was needed before) waveforms. + * @param[in] fragment_wvfms A 2D vector containing the waveforms from the current fragment to be combined. + * @param[in] num_channels The number of channels per board. + * + * @details + * The function performs the following steps: + * 1. Checks if the `wvfms` vector is empty. If it is, it resizes it to accommodate `num_channels` channels. + * 2. Iterates over each channel from 0 to `num_channels - 1`. + * 3. For each channel, it appends the samples from `fragment_wvfms` to the corresponding channel in `wvfms`. + * + * This approach ensures that waveforms from multiple fragments are concatenated correctly, maintaining the order of samples for each channel. + * + * @pre The `fragment_wvfms` vector should contain waveforms for all channels of the board being processed. + * @pre The `wvfms` vector should be either empty or already contain waveforms for all channels of the board being processed. + */ +void sbndaq::SBNDXARAPUCADecoder::append_waveforms(std::vector>& wvfms, const std::vector>& fragment_wvfms, uint32_t num_channels) { + if (fdebug_extended_fragments) std::cout << " > SBNDXARAPUCADecoder::append_waveforms: combining waveforms from extended fragments..." << std::endl; + if (wvfms.empty()) { + if (fdebug_extended_fragments) { + std::cout << "\t\t Empty waveforms, resizing to " << num_channels << " channels." << std::endl; + } + wvfms.resize(num_channels); + if (fverbose | fdebug_extended_fragments) { + std::cout << "\t\t Waveforms size BEFORE combining = " << wvfms.size() << " x " << wvfms[0].size() << std::endl; } + } + for (uint32_t ch = 0; ch < num_channels; ch++) { + wvfms[ch].insert(wvfms[ch].end(), fragment_wvfms[ch].begin(), fragment_wvfms[ch].end()); + } + if (fverbose | fdebug_extended_fragments) { + std::cout << "\t\t Waveforms size AFTER combining = " << wvfms.size() << " x " << wvfms[0].size() << std::endl; + } +} + +// =============== Dumps the decoded waveforms =============== // + +/** + * @brief Dump products into a `raw::OpDetWaveform` object and into a debug histogram file. + * @param[in,out] prod_wvfms A reference to the vector where the produced `raw::OpDetWaveform` objects are dumped into products. + * @param[in,out] wvfms A 2D vector containing the (combined if needed) waveforms. + * @param[in,out] fragment_indices A reference to a vector keeping track of the number of fragments decoded per board. + * @param[in] board_idx The board index (position in the list of boards). + * @param[in] num_channels The number of channels per board. + * @param[in] ini_wvfm_timestamp The initial timestamp of the waveform in microseconds. + * @param[in] end_wvfm_timestamp The final timestamp of the waveform in microseconds. + * + * @details + * The function performs the following steps: + * 1. Determines the number of debug waveforms to be stored based on the configuration parameter `fstore_debug_waveforms`. + * - If `fstore_debug_waveforms` is set to -1, all channels are considered for debug storage. + * - Otherwise, it takes the minimum between `num_channels` and `fstore_debug_waveforms`. + * 2. Iterates over each channel up to the determined number of debug waveforms: + * - Calls `save_prod_wvfm` to convert and store the waveform in the products. + * - Calls `save_debug_wvfm` to save the waveform as a histogram for debugging purposes. + * 3. For any remaining channels beyond the debug limit, it only calls `save_prod_wvfm` to store the waveform in the products. + * 4. Clears the `wvfms` vector to free up memory after processing. + * + * @pre ini_wvfm_timestamp and end_wvfm_timestamp is assumed to be given in microseconds when the timing frame is not CAEN_ONLY_TIMING. + * @pre The `wvfms` vector should contain waveforms for all channels of the board specified by `board_idx`. + * @pre The `fragment_indices` vector should have been initialized with a size equal to the number of boards being processed. + * @pre The `prod_wvfms` vector should be ready to accept new `raw::OpDetWaveform` objects. + * + * @see raw::OpDetWaveform + * @see save_prod_wvfm + * @see save_debug_wvfm + */ +void sbndaq::SBNDXARAPUCADecoder::dump_waveforms(std::vector & prod_wvfms, std::vector>& wvfms, size_t first_frag_idx, size_t board_idx, uint32_t num_channels, double ini_wvfm_timestamp, double end_wvfm_timestamp) { + + // The decoded waveforms are dumped into two products: + // - A xarapucadecoder-art.root file with the OpDetWaveforms as the product of this producer for further analysis. + // - A decoder_hist.root file gathering a waveform histograms. + if (fverbose) std::cout << " > SBNDXARAPUCADecoder::decode_fragment: binary decoding complete, dumping products..." << std::endl; + + uint32_t num_debug_wvfms; + + if (fstore_debug_waveforms == -1) { + num_debug_wvfms = num_channels; + } else { + num_debug_wvfms = std::min(num_channels, fstore_debug_waveforms); + } + + uint32_t ch; + + for (ch = 0; ch < num_debug_wvfms; ch++) { + save_prod_wvfm(board_idx, ch, ini_wvfm_timestamp, wvfms, prod_wvfms); + save_debug_wvfm(board_idx, first_frag_idx, ch, ini_wvfm_timestamp, end_wvfm_timestamp, wvfms); + } - fragment_indices[board_idx]++; + for (;ch < num_channels; ch++) { + save_prod_wvfm(board_idx, ch, ini_wvfm_timestamp, wvfms, prod_wvfms); } + + wvfms.clear(); + } /** @@ -842,7 +912,7 @@ void sbndaq::SBNDXARAPUCADecoder::save_prod_wvfm(size_t board_idx, size_t ch, do * * @details * The function generates a histogram with: - * - X-axis representing time in microseconds if the timing frame is SPEC_TDC_TIMING or PTB_TIMING, otherwise in samples. + * - X-axis representing time in microseconds if the timing frame is from the FRAMESHIFT reference, otherwise in samples. * - Y-axis representing ADC counts. * The histogram is stored using ROOT's `TFileService`. * @@ -858,11 +928,9 @@ void sbndaq::SBNDXARAPUCADecoder::save_debug_wvfm(size_t board_idx, size_t fragm if (factive_timing_frame != CAEN_ONLY_TIMING) { TH1I* hist = tfs->make(hist_name.str().c_str(), hist_name.str().c_str(), wvfms[ch].size(), ini_wvfm_timestamp, end_wvfm_timestamp); hist->GetYaxis()->SetTitle("ADCs"); - if (factive_timing_frame == SPEC_TDC_TIMING) { - hist->GetXaxis()->SetTitle("Time w.r.t. ETRIG (SPEC-TDC) [us]"); - } else if (factive_timing_frame == PTB_TIMING) { - hist->GetXaxis()->SetTitle("Time w.r.t. HLT ETRIG (PTB) [us]"); - } + if (factive_timing_frame == FRAMESHIFT_TIMING) { + hist->GetXaxis()->SetTitle("Time w.r.t. Frameshift reference [us]"); + } for (size_t i = 0; i < wvfms[ch].size(); i++) { hist->SetBinContent(i+1, wvfms[ch][i]); } @@ -877,63 +945,7 @@ void sbndaq::SBNDXARAPUCADecoder::save_debug_wvfm(size_t board_idx, size_t fragm } -/** - * @brief Extract a sample from a 64-bit buffer using the specified bit positions. - * - * @param[in] buffer An unsigned 64-bit integer which represents a temporal buffer for the read words and where the samples are extracted from. - * @param[in] msb An unsigned 32-bit integer representing the most significative bit (MSB) where the readout from the buffer paramter. - * @param[in] lsb An unsigned 32-bit integer representing the less significative bit (LSB) from we end read - * - * @details The function shifts the buffer to the right by the number of positions specified by `lsb` so that the least significant bit of the - * sample aligns with bit 0. It then applies a mask to isolate the bits between `lsb` and `msb`, inclusive. - * - * @return The extracted sample as a 16-bit unsigned integer. - */ -uint16_t sbndaq::SBNDXARAPUCADecoder::get_sample(uint64_t buffer, uint32_t msb, uint32_t lsb) { - uint64_t mask = (1U << (msb - lsb + 1)) - 1; - uint64_t sample = buffer >> lsb; - return sample & mask; -} - -/** - * @brief Read a 32-bit word from the data pointer and advances the pointer. - * - * @param[in, out] data_ptr A reference to a pointer pointing to the current position in the data. - * - * @details This function retrieves a 32-bit word from the memory location pointed to by `data_ptr`. After reading, it advances `data_ptr` to - * the next 32-bit word location. - * - * @return The 32-bit word read from the location pointed to by `data_ptr`. - */ -uint32_t sbndaq::SBNDXARAPUCADecoder::read_word(const uint32_t* & data_ptr) { - uint32_t word = *data_ptr; - data_ptr += 1; - return word; -} - -/** - * @brief Generates a unique global channel identifier using the board slot and the channel number of that board. - * - * @param[in] board Index of the board in `fboard_id_list` from which to derive the board slot. - * @param[in] board_channel The specific channel number on the given board. - * - * @details This function computes a `channel_id` by combining the board slot and the specific - * channel number on that board. - * The unique identifier `channel_id` (\f$ CH_{ID} $\f) is computed as follows: - * \f[ - * CH\_{ID} = B\_{ID} \times 100 + CH\_{B} - * \f] - * - * Where: - * - \f$ B\_{ID} \f$ is the fragment ID retrieved from `fboard_id_list` based on the slot. - * - \f$ CH\_B \f$ is the channel number on that board. - * - * @return A unique identifier for the specified channel as an unsigned integer. - */ -unsigned int sbndaq::SBNDXARAPUCADecoder::get_channel_id(unsigned int board, unsigned int board_channel) { - unsigned int channel_id = fboard_id_list[board] * 100 + board_channel; - return channel_id; -} +// ==================== Auxiliary functions ==================== // /** * @brief Returns the signed difference between two timestamps (t1 - t2). @@ -942,7 +954,7 @@ unsigned int sbndaq::SBNDXARAPUCADecoder::get_channel_id(unsigned int board, uns * @return The signed difference between the two timestamps (t1 - t2) in nanoseconds. * @pre The timestamp is given in nanoseconds and the maximum value is less than 2^63 ns. */ -int64_t sbndaq::SBNDXARAPUCADecoder::signed_difference(uint64_t t1, uint64_t t2) { +int64_t sbndaq::SBNDXARAPUCADecoder::signed_difference_UTC_timestamp(uint64_t t1, uint64_t t2) { int64_t signed_diff = 0; if (t2 < t1) { @@ -974,11 +986,11 @@ uint64_t sbndaq::SBNDXARAPUCADecoder::abs_difference(uint64_t t1, uint64_t t2) { } /** -* @brief Formats a timestamp in nanoseconds into a easily readable string format. -* @param[in] timestamp The timestamp in nanoseconds to be formatted. -* @return A string representation of the timestamp in the format "(seconds)nanoseconds ns". -* @details The function divides the input timestamp by 1,000,000,000 to obtain the seconds component and uses the modulus operator to get the remaining nanoseconds. -*/ + * @brief Formats a timestamp in nanoseconds into a easily readable string format. + * @param[in] timestamp The timestamp in nanoseconds to be formatted. + * @return A string representation of the timestamp in the format "(seconds)nanoseconds ns". + * @details The function divides the input timestamp by 1,000,000,000 to obtain the seconds component and uses the modulus operator to get the remaining nanoseconds. + */ std::string sbndaq::SBNDXARAPUCADecoder::print_timestamp(uint64_t timestamp) { return "(" + std::to_string(timestamp / NANOSEC_IN_SEC) + ")" + std::to_string(timestamp % NANOSEC_IN_SEC) + " ns"; } diff --git a/sbndcode/Decoders/XARAPUCA/run_xarapucadecoder.fcl b/sbndcode/Decoders/XARAPUCA/run_xarapucadecoder.fcl index 710a35c71..8f89b66d6 100644 --- a/sbndcode/Decoders/XARAPUCA/run_xarapucadecoder.fcl +++ b/sbndcode/Decoders/XARAPUCA/run_xarapucadecoder.fcl @@ -1,6 +1,5 @@ #include "xarapucadecoder.fcl" -#include "spectdc_decoder_sbnd.fcl" -#include "SBNDPTBDecoderDefaults.fcl" +#include "frameshift_sbnd_data.fcl" process_name: XARAPUCADecoder @@ -24,8 +23,7 @@ outputs: { physics: { producers: { - tdcdecoder: @local::spec_tdc_decoder_sbnd - ptbdecoder: @local::SBNDPTBDecoderDefaults + frameshift: @local::frameshift_data xarapucadecoder: @local::xarapucadecoder } @@ -34,7 +32,7 @@ physics: { analyzers:{} - reco: [tdcdecoder, ptbdecoder, xarapucadecoder] + reco: [frameshift, xarapucadecoder] # Defines the output stream, there could be more than one if using filters. diff --git a/sbndcode/Decoders/XARAPUCA/xarapucadecoder.fcl b/sbndcode/Decoders/XARAPUCA/xarapucadecoder.fcl index 3456450f9..0351ba27f 100644 --- a/sbndcode/Decoders/XARAPUCA/xarapucadecoder.fcl +++ b/sbndcode/Decoders/XARAPUCA/xarapucadecoder.fcl @@ -14,32 +14,26 @@ xarapucadecoder: fragment_id_offset: 41216 # Offset to be applied to the fragment IDs (0xA100). board_id_list: [7, 13, 16, 19] # Slots where the boards are installed. # - Timing configuration. - timing_priority: 0 # 0 SPEC-TDC ETRIG, 1 PTB HLT (1-5, 14,15), 2 CAEN-only. ns_per_sample: 16 # Nanoseconds per sample. - # - RAW timestamp configuration: - raw_event_header_product_name: "daq" # Label identifying the module that created the artdaq::detail::RawEventHeader data product. - raw_trig_max_diff: 3000000 # Maximum difference (in ns) allowed between trigger timestamps and raw timestamp. - raw_timestamp_correction: 367000 # Correction (in ns) to be applied to the raw timestamp. - # - SPEC-TDC access configuration - spectdc_product_name: "tdcdecoder" # Name for getting SPEC-TDC Decoder products (if any). - spectdc_ftrig_ch: 3 # Channel assigned to flash triggers. - spectdc_etrig_ch: 4 # Channel assigned to ETT (Event Trigger Timestamp) triggers. - # - PTB access configuration - ptb_product_name: "ptbdecoder" # Name for getting PTB Decoder products (if any). - allowed_hl_triggers: [1, 2, 3, 4, 5, 14, 15] # List of allowed HLT trigger types (1-5, 14, 15). + # - FrameShift access configuration + frameshift_product_name: "frameshift" # Name for getting FrameShift module products (if any). # OUTPUT configuration: # - Output data product instance name. waveforms_instance_name: "XARAPUCAChannels" # Name for the instance product containing the raw decoded waveforms. timing_ref_instance_name: "" # Name for the instance product containing the timing reference information. store_debug_waveforms: 0 # Number of waveforms to store (0: none, -1: all, n: first n waveforms). + # Combination of extended fragments. + combine_ext_frag: true # (De)activates the combination of extended fragments into a single raw::OpDetWaveform object. + allowed_jittering: 64 # Allowed jittering (in ns) between fragments to be combined. # - Debug options. - debug_tdc_handle: false # (De)activates SPEC-TDC art::Handle information printing. - debug_ptb_handle: false # (De)activates PTB art::Handle information printing. + debug_frameshift_handle: false # (De)activates FrameShift art::Handle information printing. debug_fragments_handle: false # (De)activates V1740B CAEN fragments art::Handle information printing. debug_timing: false # (De)activates timing data printing. debug_buffer: false # (De)activates buffer status printing. debug_waveforms: false # (De)activates waveforms decoding printing. + debug_extended_fragments: false # (De)activates extended fragments information printing. + debug_jittering: false # (De)activates trigger jittering information printing. # - Verbose option. verbose: false # (De)activates verbosity. } diff --git a/sbndcode/Decoders/run_decoders_job.fcl b/sbndcode/Decoders/run_decoders_job.fcl index adb9f661d..a4fa855f5 100644 --- a/sbndcode/Decoders/run_decoders_job.fcl +++ b/sbndcode/Decoders/run_decoders_job.fcl @@ -37,6 +37,7 @@ #include "crt_decoder_sbnd.fcl" #include "crt_channel_map_service_sbnd.fcl" +#include "frameshift_sbnd_data.fcl" process_name: DECODE @@ -61,6 +62,7 @@ physics: xarapucadecoder: @local::xarapucadecoder ptbdecoder: @local::SBNDPTBDecoderDefaults crtdecoder: @local::crt_decoder + frameshift: @local::frameshift_data } filters: @@ -72,8 +74,8 @@ physics: analyzers: {} - dropchoppy: [ daq, ntbdecoder, crtdecoder, ptbdecoder, tdcdecoder, pmtdecoder, xarapucadecoder, filterchoppy ] - savechoppy: [ daq, ntbdecoder, crtdecoder, ptbdecoder, tdcdecoder, pmtdecoder, xarapucadecoder, "!filterchoppy" ] + dropchoppy: [ daq, ntbdecoder, crtdecoder, ptbdecoder, tdcdecoder, frameshift, pmtdecoder, xarapucadecoder, filterchoppy ] + savechoppy: [ daq, ntbdecoder, crtdecoder, ptbdecoder, tdcdecoder, frameshift, pmtdecoder, xarapucadecoder, "!filterchoppy" ] trigger_paths: [ dropchoppy, savechoppy ] stream1: [ out1 ] diff --git a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl index 13452e2ff..32c9d7be8 100755 --- a/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl +++ b/sbndcode/JobConfigurations/standard/reco/reco2_data.fcl @@ -1,7 +1,6 @@ #include "opt0finder_sbnd_data.fcl" #include "sbnd_tpcpmt3dbarycentermatching_config.fcl" #include "sbnd_lightpropagationcorrection_config.fcl" -#include "frameshift_sbnd_data.fcl" #include "standard_reco2_sbnd.fcl" services: @@ -33,15 +32,13 @@ physics.producers: ### SCE-aware shower reco for data pandoraSCEShower: @local::sbnd_sce_incremental_pandoraModularShowerCreationData pandoraSCEShowerSBN: @local::sbnd_sce_sbn_pandoraModularShowerCreationData - - frameshift: @local::frameshift_data } physics.reco2: [ pandora, pandoraTrack, blipPandoraCopy, blipPandoraTrackCopy, pandoraShower, pandoraShowerSBN, pandoraCaloData, pandoraPidData, pandoraLikePidData, cvn, opt0finder, crtveto, crtspacepointmatching, crttrackmatching, tpcpmtbarycentermatching, pandoraSCE, pandoraSCETrack, pandoraSCEShower, pandoraSCEShowerSBN, pandoraSCECaloData, pandoraSCEPidData, pandoraSCELikePidData, cvnSCE, opt0finderSCE, tpcpmtbarycentermatchingSCE, crtspacepointmatchingSCE, crttrackmatchingSCE, - caloskimCalorimetry, blipreco, lightpropagationcorrectionSCE, frameshift] + caloskimCalorimetry, blipreco, lightpropagationcorrectionSCE] physics.analyzers.caloskim.G4producer: "" diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc index 4905589f3..99724616e 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.cc @@ -13,7 +13,6 @@ sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet fOpFlashLabel_tpc1 ( p.get("OpFlashLabel_tpc1") ), fSpacePointLabel( p.get("SpacePointLabel") ), fOpHitsModuleLabel( p.get("OpHitsModuleLabel") ), - fSPECTDCLabel( p.get("SPECTDCLabel") ), fFlashMatchingTool( p.get("FlashMatchingTool") ), fSaveCorrectionTree( p.get("SaveCorrectionTree") ), fSpeedOfLight( p.get("SpeedOfLight") ), @@ -25,12 +24,13 @@ sbnd::LightPropagationCorrection::LightPropagationCorrection(fhicl::ParameterSet // // More initializers here. { + fNOpChannels = fWireReadout.NOpChannels(); // Initialize the TimeCorrectionVector PerChannel - for(size_t i = 0; i < fWireReadout.NOpChannels(); ++i) { + for(size_t i = 0; i GetPropagationTimeCorrectionPerChannel(); - // Get the SPECTDC product required to go to the RWM reference frame - - art::Handle> tdcHandle; - e.getByLabel(fSPECTDCLabel, tdcHandle); - if (!tdcHandle.isValid() || tdcHandle->size() == 0){ - std::cout << "No SPECTDC products found. Skip this event." << std::endl; - ResetSliceInfo(); - continue; - } - else{ - const std::vector tdc_v(*tdcHandle); - for (size_t i=0; i> flashFM; if(fFlashMatchingTool == "OpT0Finder" ){ @@ -292,10 +273,10 @@ void sbnd::LightPropagationCorrection::produce(art::Event & e) 100., -1., Ycenter, Ywidth, Zcenter, Zwidth); newFlashTime = flasht0; sbn::CorrectedOpFlashTiming correctedOpFlashTiming; - correctedOpFlashTiming.OpFlashT0 = originalFlashTime + fEventTriggerTime/1000 - fRWMTime/1000; + correctedOpFlashTiming.OpFlashT0 = originalFlashTime; correctedOpFlashTiming.NuToFLight = (Zcenter/fSpeedOfLight)/1000; correctedOpFlashTiming.NuToFCharge = (fRecoVz/fSpeedOfLight)/1000; - correctedOpFlashTiming.OpFlashT0Corrected = newFlashTime + fEventTriggerTime/1000 - fRWMTime/1000; + correctedOpFlashTiming.OpFlashT0Corrected = newFlashTime; correctedOpFlashTimes->emplace_back(std::move(correctedOpFlashTiming)); } @@ -322,8 +303,6 @@ void sbnd::LightPropagationCorrection::beginJob() fTree->Branch("eventID", &fEvent, "eventID/i"); fTree->Branch("runID", &fRun, "runID/i"); fTree->Branch("subrunID", &fSubrun, "subrunID/i"); - fTree->Branch("RWMTime", &fRWMTime); - fTree->Branch("EventTriggerTime", &fEventTriggerTime); fTree->Branch("NuScore", &fNuScore); fTree->Branch("FMScore", &fFMScore); fTree->Branch("OpFlashTimeOld", &fOpFlashTimeOld); @@ -357,8 +336,6 @@ void sbnd::LightPropagationCorrection::ResetEventVars() fRun = 0; fSubrun = 0; _fNuScore = 0.0; - fRWMTime=-99999.; - fEventTriggerTime=-99999.; fNuScore.clear(); fFMScore.clear(); fOpFlashTimeOld.clear(); @@ -418,7 +395,7 @@ void sbnd::LightPropagationCorrection::ResetSliceInfo() fSpacePointY.clear(); fSpacePointZ.clear(); fSpacePointIntegral.clear(); - fTimeCorrectionPerChannel.resize(312, 0.0); // Reset the time correction vector for each channel + fTimeCorrectionPerChannel.resize(fNOpChannels, 0.0); // Reset the time correction vector for each channel fChargeBarycenterX.assign(2, 0.0); fChargeBarycenterY.assign(2, 0.0); fChargeBarycenterZ.assign(2, 0.0); @@ -442,16 +419,16 @@ void sbnd::LightPropagationCorrection::GetPropagationTimeCorrectionPerChannel() double dy = fSpacePointY[sp] - _opDetY; double dz = fSpacePointZ[sp] - _opDetZ; double distanceToOpDet = std::sqrt(dx*dx + dy*dy + dz*dz); - //double spToCathode = abs(fSpacePointX[sp]); // Distance from space point to cathode in mm - //double cathodeToOpDet = std::sqrt(_opDetX*_opDetX + dy*dy + dz*dz); // Distance from cathode to OpDet in mm - //float lightPropTimeVIS = spToCathode/fVGroupVUV + cathodeToOpDet/fVGroupVIS; // Speed - double cathodeToOpDet = std::sqrt(_opDetX*_opDetX + (dy/2)*(dy/2) + (dz/2)*(dz/2)); // Distance from cathode to OpDet in mm double spToCathode = std::sqrt( fSpacePointX[sp]*fSpacePointX[sp] + (dy/2)*(dy/2) + (dz/2)*(dz/2)); // Distance from space point to cathode in mm float lightPropTimeVIS = spToCathode/fVGroupVUV + cathodeToOpDet/fVGroupVIS; // Speed float lightPropTimeVUV = distanceToOpDet / fVGroupVUV; // Speed of light in mm/ns for VUV - float lightPropTime = std::min(lightPropTimeVIS, lightPropTimeVUV); + float lightPropTime = 0; + if(fPDSMap.pdType(opdet)=="pmt_coated" || fPDSMap.pdType(opdet)=="xarapuca_vuv") + lightPropTime = std::min(lightPropTimeVIS, lightPropTimeVUV); + else if(fPDSMap.pdType(opdet)=="pmt_uncoated" || fPDSMap.pdType(opdet)=="xarapuca_vis") + lightPropTime = lightPropTimeVIS; float partPropTime = std::sqrt((fSpacePointX[sp]-fRecoVx)*(fSpacePointX[sp]-fRecoVx) + (fSpacePointY[sp]-fRecoVy)*(fSpacePointY[sp]-fRecoVy) + (fSpacePointZ[sp]-fRecoVz)*(fSpacePointZ[sp]-fRecoVz))/fSpeedOfLight; float PropTime = lightPropTime + partPropTime; if(PropTime < minPropTime) minPropTime = PropTime; @@ -543,4 +520,4 @@ ::lightana::LiteOpHitArray_t sbnd::LightPropagationCorrection::GetAssociatedLite } return flash_hits_v; - } \ No newline at end of file + } diff --git a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh index e5427282a..d0ae118fd 100644 --- a/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh +++ b/sbndcode/LightPropagationCorrection/LightPropagationCorrection_module.hh @@ -145,7 +145,7 @@ private: geo::WireReadoutGeom const& fWireReadout = art::ServiceHandle()->Get(); - + opdet::sbndPDMapAlg fPDSMap; //Flash finder manager ::lightana::FlashFinderManager _mgr; ::lightana::FlashFinderManager _mgr_tpc0; @@ -171,7 +171,6 @@ private: std::string fSpacePointLabel; std::string fOpHitsModuleLabel; std::string fOpFlashNewLabel; - std::string fSPECTDCLabel; std::string fFlashMatchingTool; bool fSaveCorrectionTree; @@ -215,10 +214,9 @@ private: int fEvent; int fRun; int fSubrun; + size_t fNOpChannels; double _fNuScore; double _fFMScore; - double fEventTriggerTime=-999999.; - double fRWMTime=-999999.; std::vector fNuScore; std::vector fFMScore; std::vector fOpFlashTimeOld; diff --git a/sbndcode/Timing/CMakeLists.txt b/sbndcode/Timing/CMakeLists.txt index 1ae0bbd44..9db5843cf 100644 --- a/sbndcode/Timing/CMakeLists.txt +++ b/sbndcode/Timing/CMakeLists.txt @@ -1,9 +1,14 @@ -art_dictionary(DICTIONARY_LIBRARIES -lardataobj::RawData) +art_make_library( + SOURCE TimingUtils.cc + LIBRARIES art::Persistency_Common + messagefacility::MF_MessageLogger + ) + +art_dictionary(DICTIONARY_LIBRARIES lardataobj::RawData) install_headers() -# install_fhicl() -# install_source() +#install_fhicl() +install_source() add_subdirectory(FrameShift) add_subdirectory(WaveformAlignment) diff --git a/sbndcode/Timing/FrameShift/FrameShift_module.cc b/sbndcode/Timing/FrameShift/FrameShift_module.cc index c675eb5a8..62e9d5a27 100644 --- a/sbndcode/Timing/FrameShift/FrameShift_module.cc +++ b/sbndcode/Timing/FrameShift/FrameShift_module.cc @@ -36,12 +36,12 @@ #include -namespace sbnd { +namespace sbnd::timing { class FrameShift; } -class sbnd::FrameShift : public art::EDProducer { +class sbnd::timing::FrameShift : public art::EDProducer { public: explicit FrameShift(fhicl::ParameterSet const& p); // The compiler-generated destructor is fine for non-base @@ -59,708 +59,676 @@ class sbnd::FrameShift : public art::EDProducer { // Selected optional functions. void beginJob() override; void ResetEventVars(); - double SubtractUTCTimestmap(const uint64_t& ts1, const uint64_t& ts2); -private: + uint64_t FindClosest(const std::vector ×tamps, const uint64_t &reference); + std::string PrintFormatTimestamp(const uint64_t ×tamp); - // Declare member data here. - //DAQ Header - std::string fDAQHeaderModuleLabel; - std::string fDAQHeaderInstanceLabel; - uint64_t fRawTSCorrection = 0; //Correction for the Event timestamp to account for NTB - uint64_t fMaxAllowedRefTimeDiff; - uint64_t _raw_ts = 0; //ns - - //Timing Reference - art::InputTag fTimingRefPmtLabel; - art::InputTag fTimingRefCrtLabel; + void GetRawTimestamp(const art::Event &e); + void GetTDCTimestamps(const art::Event &e); + void GetPTBTimestamps(const art::Event &e); + void FindETRIGs(); - uint16_t _pmt_timing_type; // e.g. SPECTDC = 0; PTB HLT = 1; Do nothing = 2 - uint16_t _pmt_timing_ch; - uint16_t _crt_timing_type; - uint16_t _crt_timing_ch; + uint64_t DecideGlobalEtrigTimestamp(); + void DecideRelevantTDCTimestamps(const uint64_t &global_etrig_ts); + void DecideRelevantPTBTimestamps(const uint64_t &global_etrig_ts); - //PTB - art::InputTag fPtbDecodeLabel; - - std::vector fPtbEtrigHlts; - - std::vector fBeamEtrigHlt; - std::vector fOffbeamEtrigHlt; - std::vector fXmuonEtrigHlt; - - int fBeamCrtT1Hlt; - int fOffbeamCrtT1Hlt; - - int fBeamGateHlt; - int fOffbeamGateHlt; - - std::vector _ptb_hlt_trigger; - std::vector _ptb_hlt_timestamp; - std::vector _ptb_hlt_unmask_timestamp; - std::vector _ptb_hlt_trunmask; - - std::vector _ptb_hlt_etrig; - std::vector _ptb_hlt_etrig_ts; - - bool _isBeam; - bool _isOffbeam; - bool _isXmuon; - int _hlt_etrig; - uint64_t _hlt_etrig_ts; - int _hlt_gate; - uint64_t _hlt_gate_ts; - int _hlt_crtt1; - uint64_t _hlt_crtt1_ts; +private: + + // FCL Parameters + std::string fDAQHeaderModuleLabel; // Instance label for DAQ header product + std::string fDAQHeaderInstanceLabel; // Module label for DAQ header product + uint64_t fRawTSCorrection; // Correction to the DAQ header timestamp to account for NTB offset + art::InputTag fPtbDecodeLabel; // Module label for decoded PTB data + std::vector fPtbEtrigHlts; // All allowed PTB HLT codes for event triggers + std::vector fBeamEtrigHlt; // PTB HLT codes associated with beam trigger + std::vector fOffbeamEtrigHlt; // PTB HLT codes associated with offbeam trigger + std::vector fXmuonEtrigHlt; // PTB HLT codes associated with crossing muon trigger + int fBeamCrtT1Hlt; // PTB HLT code for beam CRT T1 reset + int fOffbeamCrtT1Hlt; // PTB HLT code for offbeam CRT T1 reset + int fBeamGateHlt; // PTB HLT code for beam gate opening + int fOffbeamGateHlt; // PTB HLT code for offbeam gate opening + art::InputTag fTdcDecodeLabel; // Module label for decoded TDC data + uint64_t fShiftData2MC; // Value to shift Data to MC -- so that data agree with MC [ns] + // TODO: Derive this value and verify if it is consistent across pmt/crt + // TODO: Get this value from database instead of fhicl parameter + uint64_t fShiftTdcRwm2PtbGate; // Value to move Tdc RWM (ch2) frame to agree with HLT Gate Frame + // TODO: Get this value from database instead of fhicl parameter + uint64_t fShiftTdcEtrig2PtbEtrig; // Value to move Tdc Etrig (ch4) to agree with HLT ETRIG + // TODO: Get from database + uint64_t fShiftTdcCrtt12PtbCrtt1; // Value to move Tdc Crt T1 Reset (ch0) to agree with HLT CRT T1 + // TODO: Get from database + bool fMakeTree; // Whether to produce a TTree in the hist file + bool fDebugDAQHeader; // Whether to print debug statements relevant to DAQ header + bool fDebugPtb; // Whether to print debug statements relevant to PTB + bool fDebugTdc; // Whether to print debug statements relevant to TDC + bool fDebugFrame; // Whether to print debug statements relevant to frame shifts + + // Global Variables, set in processing + int _run, _subrun, _event; // Stores the unique run, subrun, event number combination for this event + + uint64_t _raw_ts; // Stores DAQ header timestamp + std::vector _ptb_hlt_trigger; // Stores full trigger word from PTB HLT (indexed per HLT) + std::vector _ptb_hlt_timestamp; // Stores timestamp associated with PTB HLT (indexed per HLT) + std::vector _ptb_hlt_unmask_timestamp; // Stores timestamp associated with PTB HLT for each element (indexed per HLT element once unmasked) + std::vector _ptb_hlt_trunmask; // Stores PTB HLT code for each element (indexed per HLT element once unmasked) + std::vector _ptb_hlt_etrig; // Stores all PTB HLTs that are an 'allowed type' for this stream + std::vector _ptb_hlt_etrig_ts; // Stores the associated timestamps for the above HLTs + + bool _isBeam; // Event is a beam trigger (BNB or BNBLight) + bool _isOffbeam; // Event is an offbeam trigger (OffBeam or OffBeamLight) + bool _isXmuon; // Event is a crossing muon trigger + + int _hlt_etrig; // Stores the HLT event trigger code, once we've decided which one was closest to the DAQ header timestamp + uint64_t _hlt_etrig_ts; // Stores the HLT timestamp, once we've decided which one was closest to the DAQ header timestamp + int _hlt_gate; // Stores the HLT gate opening code, once we've decided which one was closest to the DAQ header timestamp + uint64_t _hlt_gate_ts; // Stores the HLT gate opening timestamp, once we've decided which one was closest to the DAQ header timestamp + int _hlt_crtt1; // Stores the HLT CRT T1 code, once we've decided which one was closest to the DAQ header timestamp + uint64_t _hlt_crtt1_ts; // Stores the HLT CRT T1 timestamp, once we've decided which one was closest to the DAQ header timestamp - //TDC - art::InputTag fTdcDecodeLabel; - - std::vector _tdc_ch0; - std::vector _tdc_ch1; - std::vector _tdc_ch2; - std::vector _tdc_ch4; - - uint64_t _tdc_crtt1_ts; - uint64_t _tdc_bes_ts; - uint64_t _tdc_rwm_ts; - uint64_t _tdc_etrig_ts; - - - //Frame Shift - uint16_t _global_timing_type; // e.g. SPECTDC = 0; PTB HLT = 1; - uint64_t _global_frame; - - double _frame_tdc_crtt1; //ns - double _frame_tdc_bes; //ns - double _frame_tdc_rwm; //ns - double _frame_hlt_crtt1; //ns - double _frame_hlt_gate; //ns - - //Value to apply at CAF depending on which stream - double _frame_apply_at_caf; - - //Value to shift Data to MC -- so that data agree with MC [ns] - //TODO: Derive this value and verify if it is consistent across pmt/crt - //TODO: Get this value from database instead of fhicl parameter - double _frame_data2mc; //ns - - //Value to move RWM frame to agree with HLT Gate Frame - //This is derived by subtracting: TDC RWM - HLT Gate. - //Using the MC2025B dataset, this distribution has a mean of 1738 ns and std of 9 ns. - //TODO: Get this value from database instead of fhicl parameter - double fShiftRWM2Gate; //ns - - //Debug - bool fDebugDAQHeader; - bool fDebugTimingRef; - bool fDebugPtb; - bool fDebugTdc; - bool fDebugFrame; - - //---TREE PARAMETERS + std::vector _tdc_ch0; // Stores all the timestamps recorded in the TDC channel 0 (CRT T1) + std::vector _tdc_ch1; // Stores all the timestamps recorded in the TDC channel 1 (BES) + std::vector _tdc_ch2; // Stores all the timestamps recorded in the TDC channel 2 (RWM) + //std::vector _tdc_ch3; // Stores all the timestamps recorded in the TDC channel 3 (Flash trigger) + std::vector _tdc_ch4; // Stores all the timestamps recorded in the TDC channel 4 (Event trigger) + + uint64_t _tdc_crtt1_ts; // Stores the TDC CRT T1 timestamp, once we've decided which one was closest to the DAQ header timestamp + uint64_t _tdc_bes_ts; // Stores the TDC BES timestamp, once we've decided which one was closest to the DAQ header timestamp + uint64_t _tdc_rwm_ts; // Stores the TDC RWM timestamp, once we've decided which one was closest to the DAQ header timestamp + uint64_t _tdc_etrig_ts; // Stores the TDC event trigger timestamp, once we've decided which one was closest to the DAQ header timestamp + + uint64_t _frame_crtt1; // Stores the shift required to move from the PPS frame to the CRT T1 frame + uint16_t _timing_type_crtt1; // Stores the type (TDC/PTB) used to produce the above shift + uint16_t _timing_channel_crtt1; // Stores the channel (TDC) / code (PTB) used to produce the above shift + uint64_t _frame_gate; // Stores the shift required to move from the PPS frame to the gate opening frame + uint16_t _timing_type_gate; // Stores the type (TDC/PTB) used to produce the above shift + uint16_t _timing_channel_gate; // Stores the channel (TDC) / code (PTB) used to produce the above shift + uint64_t _frame_etrig; // Stores the shift required to move from the PPS frame to the event trigger frame + uint16_t _timing_type_etrig; // Stores the type (TDC/PTB) used to produce the above shift + uint16_t _timing_channel_etrig; // Stores the channel (TDC) / code (PTB) used to produce the above shift + + uint64_t _frame_default; // Stores the shift required to move from the PPS frame to the default frame for this stream + uint16_t _timing_type_default; // Stores the type (TDC/PTB) used to produce the above shift + uint16_t _timing_channel_default; // Stores the channel (TDC) / code (PTB) used to produce the above shift + + // Tree production TTree *fTree; art::ServiceHandle tfs; - int _run, _subrun, _event; }; -sbnd::FrameShift::FrameShift(fhicl::ParameterSet const& p) - : EDProducer{p} // , - // More initializers here. +sbnd::timing::FrameShift::FrameShift(fhicl::ParameterSet const& p) + : EDProducer{p} { fDAQHeaderInstanceLabel = p.get("DAQHeaderInstanceLabel"); fDAQHeaderModuleLabel = p.get("DAQHeaderModuleLabel"); fRawTSCorrection = p.get("RawTSCorrection"); - fMaxAllowedRefTimeDiff = p.get("MaxAllowedRefTimeDiff"); - - fTimingRefPmtLabel = p.get("TimingRefPmtLabel"); - fTimingRefCrtLabel = p.get("TimingRefCrtLabel"); - - fTdcDecodeLabel = p.get("TdcDecodeLabel"); fPtbDecodeLabel = p.get("PtbDecodeLabel"); - fPtbEtrigHlts = p.get>("PtbEtrigHlts"); - fBeamEtrigHlt = p.get>("BeamEtrigHlt"); fOffbeamEtrigHlt = p.get>("OffbeamEtrigHlt"); fXmuonEtrigHlt = p.get>("XmuonEtrigHlt"); - fBeamCrtT1Hlt = p.get("BeamCrtT1Hlt"); fOffbeamCrtT1Hlt = p.get("OffbeamCrtT1Hlt"); - fBeamGateHlt = p.get("BeamGateHlt"); fOffbeamGateHlt = p.get("OffbeamGateHlt"); - + fTdcDecodeLabel = p.get("TdcDecodeLabel"); + fShiftData2MC = p.get("ShiftData2MC"); + fShiftTdcRwm2PtbGate = p.get("ShiftTdcRwm2PtbGate"); + fShiftTdcEtrig2PtbEtrig = p.get("ShiftTdcEtrig2PtbEtrig"); + fShiftTdcCrtt12PtbCrtt1 = p.get("ShiftTdcCrtt12PtbCrtt1"); + fMakeTree = p.get("MakeTree", false); fDebugDAQHeader = p.get("DebugDAQHeader", false); - fDebugTimingRef = p.get("DebugTimingRef", false); fDebugPtb = p.get("DebugPtb", false); fDebugTdc = p.get("DebugTdc", false); fDebugFrame = p.get("DebugFrame", false); - //TODO: Get from database instead of fhicl parameters - _frame_data2mc = p.get("ShiftData2MC"); - - //TODO: Get from database instead of fhicl parameters - fShiftRWM2Gate = p.get("ShiftRWM2Gate"); - - produces< sbnd::timing::FrameShiftInfo >(); - produces< sbnd::timing::TimingInfo >(); + produces(); + produces(); } -void sbnd::FrameShift::produce(art::Event& e) +void sbnd::timing::FrameShift::produce(art::Event& e) { - ResetEventVars(); _run = e.id().run(); _subrun = e.id().subRun(); _event = e.id().event(); - if (fDebugTdc | fDebugPtb | fDebugTimingRef | fDebugFrame) - std::cout <<"#----------RUN " << _run << " SUBRUN " << _subrun << " EVENT " << _event <<"----------#\n"; - - //---------------------------TimingRef-----------------------------// - art::Handle timingRefPmtHandle; - e.getByLabel(fTimingRefPmtLabel, timingRefPmtHandle); - - if (!timingRefPmtHandle.isValid()){ - throw cet::exception("FrameShift") << "No raw::TimingReferenceInfo found w/ tag " << fTimingRefPmtLabel << ". Check data quality!"; - } - else{ - raw::TimingReferenceInfo const& pmt_timing(*timingRefPmtHandle); - _pmt_timing_type = pmt_timing.timingType; - _pmt_timing_ch = pmt_timing.timingChannel; - } - timingRefPmtHandle.removeProduct(); - - art::Handle timingRefCrtHandle; - e.getByLabel(fTimingRefCrtLabel, timingRefCrtHandle); - - if (!timingRefCrtHandle.isValid()){ - throw cet::exception("FrameShift") << "No raw::TimingReferenceInfo found w/ tag " << fTimingRefCrtLabel << ". Check data quality!"; - } - else{ - raw::TimingReferenceInfo const& crt_timing(*timingRefCrtHandle); - _crt_timing_type = crt_timing.timingType; - _crt_timing_ch = crt_timing.timingChannel; - } - timingRefCrtHandle.removeProduct(); - - //TODO: Add XA check! - - //Check what frame it is decoded to? - if (_pmt_timing_type != _crt_timing_type){ - throw cet::exception("FrameShift") << "Timing Reference for PMT and CRT are not the same! PMT type = " - << _pmt_timing_type << ", CRT type = " << _crt_timing_type; - }else{ - _global_timing_type = _pmt_timing_type; - } - - if (fDebugTimingRef){ - std::cout << "----------------------------------------------------" << std::endl; - std::cout << "Timing Reference For Decoding Pmt" << std::endl; - std::cout << " Type = " << _pmt_timing_type << " (SPECTDC = 0; PTB HLT = 1; Nothing = 2)." << std::endl; - std::cout << " Channel = " << _pmt_timing_ch << std::endl; - std::cout << "----------------------------------------------------" << std::endl; - std::cout << "Timing Reference For Decoding Crt" << std::endl; - std::cout << " Type = " << _crt_timing_type << " (SPECTDC = 0; PTB HLT = 1; Nothing = 2)." << std::endl; - std::cout << " Channel = " << _crt_timing_ch << std::endl; - std::cout << "----------------------------------------------------" << std::endl; - std::cout << "Global Timing Type = " << _global_timing_type << " (SPECTDC = 0; PTB HLT = 1; Nothing = 2)." << std::endl; - std::cout << "----------------------------------------------------" << std::endl; - } - - //----------------------------DAQ Header-----------------------------// + if(fDebugTdc | fDebugPtb | fDebugFrame) + std::cout <<"#----------RUN " << _run << " SUBRUN " << _subrun << " EVENT " << _event <<"----------#\n"; + + GetRawTimestamp(e); + GetPTBTimestamps(e); + GetTDCTimestamps(e); + FindETRIGs(); + + uint64_t global_etrig_ts = DecideGlobalEtrigTimestamp(); + DecideRelevantTDCTimestamps(global_etrig_ts); + DecideRelevantPTBTimestamps(global_etrig_ts); + + if(_isBeam || _isOffbeam) + { + //Frame CRT T1 + if(_tdc_crtt1_ts != kInvalidTimestamp) + { + _frame_crtt1 = _tdc_crtt1_ts - fShiftTdcCrtt12PtbCrtt1; + _timing_type_crtt1 = kSPECTDCType; + _timing_channel_crtt1 = 0; + } + else if(_hlt_crtt1_ts != kInvalidTimestamp) + { + _frame_crtt1 = _hlt_crtt1_ts; + _timing_type_crtt1 = kPTBHLTType; + _timing_channel_crtt1 = _hlt_crtt1; + } + else + { + _frame_crtt1 = 0; + _timing_type_crtt1 = kNoShiftType; + _timing_channel_crtt1 = kInvalidChannel; + } + + //Frame Beam Gate + if(_isBeam && _tdc_rwm_ts != kInvalidTimestamp) + { + _frame_gate = _tdc_rwm_ts - fShiftTdcRwm2PtbGate; //TODO: Add shift from Data to MD: + fShiftData2MC + _timing_type_gate = kSPECTDCType; + _timing_channel_gate = 2; + } + else if(_hlt_gate_ts != kInvalidTimestamp) + { + _frame_gate = _hlt_gate_ts; + _timing_type_gate = kPTBHLTType; + _timing_channel_gate = _hlt_gate; + } + else + { + _frame_gate = 0; + _timing_type_gate = kNoShiftType; + _timing_channel_gate = kInvalidChannel; + } + } + + //Frame ETRIG + if(_tdc_etrig_ts != kInvalidTimestamp) + { + _frame_etrig = _tdc_etrig_ts - fShiftTdcEtrig2PtbEtrig; + _timing_type_etrig = kSPECTDCType; + _timing_channel_etrig = 4; + } + else if(_hlt_etrig_ts != kInvalidTimestamp) + { + _frame_etrig = _hlt_etrig_ts; + _timing_type_etrig = kPTBHLTType; + _timing_channel_etrig = _hlt_etrig; + } + else + { + _frame_etrig = 0; + _timing_type_etrig = kNoShiftType; + _timing_channel_etrig = kInvalidChannel; + } + + //-----------------------Pick default frame-----------------------// + // The follow picks which frame to apply at downstream stage and store it as frame_default, based on the stream + // + // 1. Beam Stream: recontruct the beam spill + porch relative to RWM record in TDC since it has better resolution. + // Then, apply a constant to shift the RWM fram + // i.e. frame_tdc_rwm + shift rwm to beam gate HLT 26 + // + // 2. Offbeam Stream: reconstruct the porch relative to the beam gate opening frame + // i.e. frame_hlt_gate equiv + // + // 3. Xmuon Stream: shift to etrig + + if(_isBeam || _isOffbeam) + { + //Pick default stream -- beam gate + _frame_default = _frame_gate; + _timing_type_default = _timing_type_gate; + _timing_channel_default = _timing_channel_gate; + } + else if(_isXmuon) + { + //Pick default stream -- ETRIG + _frame_default = _frame_etrig; + _timing_type_default = _timing_type_etrig; + _timing_channel_default = _timing_channel_etrig; + } + + if(fDebugFrame) + { + std::cout << "--------------------------------------" << std::endl; + std::cout << "Frame Shift Results:" << std::endl; + std::cout << "Frame CRT T1 type " << _timing_type_crtt1 <<", channel " << _timing_channel_crtt1 << ": " + << PrintFormatTimestamp(_frame_crtt1) << std::endl; + std::cout << "Frame Beam Gate type " << _timing_type_gate <<", channel " << _timing_channel_gate << ": " + << PrintFormatTimestamp(_frame_gate) << std::endl; + std::cout << "Frame ETRIG type " << _timing_type_etrig <<", channel " << _timing_channel_etrig << ": " + << PrintFormatTimestamp(_frame_etrig) << std::endl; + std::cout << "Default Frame type " << _timing_type_default <<", channel " << _timing_channel_default << ": " + << PrintFormatTimestamp(_frame_default) << std::endl; + std::cout << "--------------------------------------" << std::endl; + } + + //Put product in event + std::unique_ptr newFrameShiftInfo(new FrameShiftInfo(_frame_crtt1, _timing_type_crtt1, _timing_channel_crtt1, + _frame_gate, _timing_type_gate, _timing_channel_gate, + _frame_etrig, _timing_type_etrig,_timing_channel_etrig, + _frame_default, _timing_type_default, _timing_channel_default)); + + std::unique_ptr newTimingInfo(new TimingInfo(_raw_ts, _tdc_crtt1_ts, _tdc_bes_ts, _tdc_rwm_ts, _tdc_etrig_ts, + _hlt_crtt1_ts, _hlt_etrig_ts, _hlt_gate_ts)); + + e.put(std::move(newTimingInfo)); + e.put(std::move(newFrameShiftInfo)); + + //Fill the tree + if(fMakeTree) + fTree->Fill(); +} + +uint64_t sbnd::timing::FrameShift::FindClosest(const std::vector ×tamps, const uint64_t &reference) +{ + uint64_t min_diff = kInvalidTimestamp; + uint64_t closest = kInvalidTimestamp; + + for(const uint64_t ×tamp : timestamps) + { + uint64_t diff = reference > timestamp ? reference - timestamp : timestamp - reference; + + if(diff < min_diff) + { + min_diff = diff; + closest = timestamp; + } + } + + return closest; +} + +std::string sbnd::timing::FrameShift::PrintFormatTimestamp(const uint64_t ×tamp) +{ + std::stringstream ss; + ss << "(s) = " << timestamp / kSecondInNanoseconds << ", (ns) = " << timestamp % kSecondInNanoseconds; + + return ss.str(); +} + +void sbnd::timing::FrameShift::GetRawTimestamp(const art::Event &e) +{ art::Handle DAQHeaderHandle; e.getByLabel(fDAQHeaderModuleLabel, fDAQHeaderInstanceLabel, DAQHeaderHandle); - if (!DAQHeaderHandle.isValid()){ + if(!DAQHeaderHandle.isValid()) throw cet::exception("FrameShift") << "No artdaq::detail::RawEventHeader found w/ tag " << fDAQHeaderModuleLabel << ". Check data quality!"; - } - else{ + else + { artdaq::RawEvent rawHeaderEvent = artdaq::RawEvent(*DAQHeaderHandle); _raw_ts = rawHeaderEvent.timestamp() - fRawTSCorrection; - } - if (fDebugDAQHeader){ - std::cout << "----------------------------------------------------" << std::endl; - std::cout << "DAQ Header Timestamp: " - << " (s) = " << _raw_ts/uint64_t(1e9) - << ", (ns) = " << _raw_ts%uint64_t(1e9) - << std::endl; - std::cout << "----------------------------------------------------" << std::endl; - } - - //---------------------------TDC-----------------------------// - art::Handle> tdcHandle; - e.getByLabel(fTdcDecodeLabel, tdcHandle); + } + + if(fDebugDAQHeader) + { + std::cout << "----------------------------------------------------" << std::endl; + std::cout << "DAQ Header Timestamp: " << PrintFormatTimestamp(_raw_ts) << std::endl; + } +} - if (!tdcHandle.isValid() || tdcHandle->size() == 0){ - mf::LogInfo("FrameShift") << "No sbnd::timing::DAQTimestamp found w/ tag " << fTdcDecodeLabel << ". Check data quality!\n"; - //throw cet::exception("FrameShift") << "No sbnd::timing::DAQTimestamp found w/ tag " << fTdcDecodeLabel << ". Check data quality!"; - } - else{ - std::vector tdc_v(*tdcHandle); - for (size_t i=0; i> tdcHandle; + e.getByLabel(fTdcDecodeLabel, tdcHandle); + + if(!tdcHandle.isValid() || tdcHandle->size() == 0) + mf::LogInfo("FrameShift") << "No DAQTimestamp found w/ tag " << fTdcDecodeLabel << ". Check data quality!\n"; + else + { + for(auto const& tdc : *tdcHandle) + { + const uint32_t ch = tdc.Channel(); + const uint64_t ts = tdc.Timestamp(); + //Also TODO: Make use of picosecond timestamps + //const uint64_t ts_ps = tdc.TimestampPs(); + + switch(ch) + { + case 0: + _tdc_ch0.push_back(ts); + break; + case 1: + _tdc_ch1.push_back(ts); + break; + case 2: + _tdc_ch2.push_back(ts); + break; + //case 3: + // _tdc_ch3.push_back(ts); + // break; + case 4: + _tdc_ch4.push_back(ts); + break; + } + } } - } - tdcHandle.removeProduct(); +} - //---------------------------PTB-----------------------------// +void sbnd::timing::FrameShift::GetPTBTimestamps(const art::Event &e) +{ art::Handle> ptbHandle; - std::vector> ptb_v; e.getByLabel(fPtbDecodeLabel, ptbHandle); - if ((!ptbHandle.isValid() || ptbHandle->size() == 0)){ - throw cet::exception("FrameShift") << "No raw::ptb::sbndptb found w/ tag " << fPtbDecodeLabel << ". Check data quality!"; - } - else{ - art::fill_ptr_vector(ptb_v, ptbHandle); - // HLT Words - unsigned nHLTs = 0; - - for(auto const& ptb : ptb_v) - nHLTs += ptb->GetNHLTriggers(); + if((!ptbHandle.isValid() || ptbHandle->size() == 0)) + throw cet::exception("FrameShift") << "No raw::ptb::sbndptb found w/ tag " << fPtbDecodeLabel << ". Check data quality!"; + else + { + // HLT Words + unsigned nHLTs = 0; + + for(auto const& ptb : *ptbHandle) + nHLTs += ptb.GetNHLTriggers(); - _ptb_hlt_trigger.resize(nHLTs); - _ptb_hlt_timestamp.resize(nHLTs); - _ptb_hlt_trunmask.resize(nHLTs); - _ptb_hlt_unmask_timestamp.resize(nHLTs); - - unsigned hlt_i = 0; //For multiple upbits in trigger words for unmasking - unsigned h_i = 0; //For trigger with bitmask - - for(auto const& ptb : ptb_v){ - for(unsigned i = 0; i < ptb->GetNHLTriggers(); ++i){ - _ptb_hlt_trigger[h_i] = ptb->GetHLTrigger(i).trigger_word; - _ptb_hlt_timestamp[h_i] = ptb->GetHLTrigger(i).timestamp * uint64_t(20); //Units can be found in the Decoder Module - h_i++; + _ptb_hlt_trigger.resize(nHLTs); + _ptb_hlt_timestamp.resize(nHLTs); + _ptb_hlt_trunmask.resize(nHLTs); + _ptb_hlt_unmask_timestamp.resize(nHLTs); + + unsigned hlt_i = 0; //For multiple upbits in trigger words for unmasking + unsigned h_i = 0; //For trigger with bitmask + + for(auto const& ptb : *ptbHandle) + { + for(unsigned i = 0; i < ptb.GetNHLTriggers(); ++i) + { + _ptb_hlt_trigger[h_i] = ptb.GetHLTrigger(i).trigger_word; + _ptb_hlt_timestamp[h_i] = ptb.GetHLTrigger(i).timestamp * uint64_t(20); //Units can be found in the Decoder Module + h_i++; - int val = ptb->GetHLTrigger(i).trigger_word; - int upBit[32]; + int val = ptb.GetHLTrigger(i).trigger_word; + int upBit[32]; - for (int u=0; u<32; u++){ //setting default values for maximum of 32 bits - upBit[u]=-1; - } + for(int u=0; u<32; u++) + upBit[u] = -1; //setting default values for maximum of 32 bits - int numOfTrig =0; - for(int b=0; b<32;b++){ - if ((val & 0x01) ==1){ - upBit[numOfTrig] = b; - numOfTrig++; - } - val = val >> 1; - } + int numOfTrig = 0; + for(int b=0; b<32; b++) + { + if((val & 0x01) == 1) + { + upBit[numOfTrig] = b; + numOfTrig++; + } + + val = val >> 1; + } - if (numOfTrig ==1){ - _ptb_hlt_unmask_timestamp[hlt_i] = _ptb_hlt_timestamp[h_i-1]; - _ptb_hlt_trunmask[hlt_i] = upBit[0]; - hlt_i++; - }//End of if statement for single upbit - else if (numOfTrig > 1){ - nHLTs += (numOfTrig -1); - _ptb_hlt_unmask_timestamp.resize(nHLTs); - _ptb_hlt_trunmask.resize(nHLTs); + if(numOfTrig == 1) + { + _ptb_hlt_unmask_timestamp[hlt_i] = _ptb_hlt_timestamp[h_i-1]; + _ptb_hlt_trunmask[hlt_i] = upBit[0]; + hlt_i++; + } + else if(numOfTrig > 1) + { + nHLTs += (numOfTrig - 1); + _ptb_hlt_unmask_timestamp.resize(nHLTs); + _ptb_hlt_trunmask.resize(nHLTs); - for (int mult =0; mult < numOfTrig; mult++){ - _ptb_hlt_trunmask[hlt_i] = upBit[mult]; - _ptb_hlt_unmask_timestamp[hlt_i] = _ptb_hlt_timestamp[h_i-1]; - hlt_i++; - } //End of loop over multiple upbits - } //End of else statement for multiple triggers - } //End of loop over nHLTriggers - } //End of loop over ptb in ptb_v - } - ptbHandle.removeProduct(); - - if (fDebugPtb){ - for (size_t i = 0; i < _ptb_hlt_unmask_timestamp.size(); i++){ - std::cout << "----------------------------------------------------" << std::endl; - std::cout << "HLT " << _ptb_hlt_trunmask[i] - << " sec (s) = " << _ptb_hlt_unmask_timestamp[i]/uint64_t(1e9) - << ", ts (ns) = " << _ptb_hlt_unmask_timestamp[i]%uint64_t(1e9) - < 1){ - //Get the one closest to the global frame - uint64_t min_diff = std::numeric_limits::max(); - for(auto ts : _tdc_ch0){ - uint64_t diff = _global_frame > ts ? _global_frame - ts : ts - _global_frame; - if(diff < min_diff) - { - min_diff = diff; - _tdc_crtt1_ts = ts; - } - } - } - - // ch1: BES - if (_tdc_ch1.size() == 1){ - _tdc_bes_ts = _tdc_ch1[0]; - } - else if(_tdc_ch1.size() > 1){ - //Get the one closest to the global frame - uint64_t min_diff = std::numeric_limits::max(); - for(auto ts : _tdc_ch1){ - uint64_t diff = _global_frame > ts ? _global_frame - ts : ts - _global_frame; - if(diff < min_diff) - { - min_diff = diff; - _tdc_bes_ts = ts; - } + for(int mult = 0; mult < numOfTrig; mult++) + { + _ptb_hlt_unmask_timestamp[hlt_i] = _ptb_hlt_timestamp[h_i-1]; + _ptb_hlt_trunmask[hlt_i] = upBit[mult]; + hlt_i++; + } + } + } //End of loop over nHLTriggers + } //End of loop over ptb in ptb_v } - } - // ch2: RWM - if (_tdc_ch2.size() == 1){ - _tdc_rwm_ts = _tdc_ch2[0]; - } - else if(_tdc_ch2.size() > 1){ - //Get the one closest to the global frame - uint64_t min_diff = std::numeric_limits::max(); - for(auto ts : _tdc_ch2){ - uint64_t diff = _global_frame > ts ? _global_frame - ts : ts - _global_frame; - if(diff < min_diff) - { - min_diff = diff; - _tdc_rwm_ts = ts; - } + if(fDebugPtb) + { + for(size_t i = 0; i < _ptb_hlt_unmask_timestamp.size(); i++) + { + std::cout << "----------------------------------------------------" << std::endl; + std::cout << "HLT " << _ptb_hlt_trunmask[i] << ": " << PrintFormatTimestamp(_ptb_hlt_unmask_timestamp[i]) << std::endl; + } } - } +} - // ch4: ETRIG - if (_tdc_ch4.size() == 0){ +void sbnd::timing::FrameShift::FindETRIGs() +{ + if(_tdc_ch4.size() == 0) throw cet::exception("FrameShift") << "No TDC ETRIG timestamps found! Check data quality!"; - } - else if (_tdc_ch4.size() == 1){ - _tdc_etrig_ts = _tdc_ch4[0]; - } - else if(_tdc_ch4.size() > 1){ - uint64_t min_diff = std::numeric_limits::max(); - for(auto ts : _tdc_ch4){ - uint64_t diff = _raw_ts > ts ? _raw_ts - ts : ts - _raw_ts; - if(diff < min_diff) - { - min_diff = diff; - _tdc_etrig_ts = ts; - } - } - } - - if (fDebugTdc){ - std::cout << "----------------------------------------------------" << std::endl; - std::cout << "TDC Channel 0 (CRTT1) Timestamp: " - << " (s) = " << _tdc_crtt1_ts/uint64_t(1e9) - << ", (ns) = " << _tdc_crtt1_ts%uint64_t(1e9) - << std::endl; - std::cout << "TDC Channel 1 (BES) Timestamp: " - << " (s) = " << _tdc_bes_ts/uint64_t(1e9) - << ", (ns) = " << _tdc_bes_ts%uint64_t(1e9) - << std::endl; - std::cout << "TDC Channel 2 (RWM) Timestamp: " - << " (s) = " << _tdc_rwm_ts/uint64_t(1e9) - << ", (ns) = " << _tdc_rwm_ts%uint64_t(1e9) - << std::endl; - std::cout << "TDC Channel 4 (ETRIG) Timestamp: " - << " (s) = " << _tdc_etrig_ts/uint64_t(1e9) - << ", (ns) = " << _tdc_etrig_ts%uint64_t(1e9) - << std::endl; - std::cout << "----------------------------------------------------" << std::endl; - } - - //---------------------------PTB Frame-----------------------------// + else + _tdc_etrig_ts = FindClosest(_tdc_ch4, _raw_ts); - //HLT ETRIG //Grab all the ETRIG HLTS -- there might be more than 1 - for (size_t i = 0; i < _ptb_hlt_unmask_timestamp.size(); i++){ - for (size_t j = 0; j < fPtbEtrigHlts.size(); j++){ - if (_ptb_hlt_trunmask[i] == fPtbEtrigHlts[j]){ - _ptb_hlt_etrig.push_back(_ptb_hlt_trunmask[i]); - _ptb_hlt_etrig_ts.push_back(_ptb_hlt_unmask_timestamp[i]); - } + for(size_t i = 0; i < _ptb_hlt_unmask_timestamp.size(); i++) + { + for(size_t j = 0; j < fPtbEtrigHlts.size(); j++) + { + if(_ptb_hlt_trunmask[i] == fPtbEtrigHlts[j]) + { + _ptb_hlt_etrig.push_back(_ptb_hlt_trunmask[i]); + _ptb_hlt_etrig_ts.push_back(_ptb_hlt_unmask_timestamp[i]); + } + } } - } - if (_ptb_hlt_etrig.size() == 0){ + if(_ptb_hlt_etrig.size() == 0) throw cet::exception("FrameShift") << "No HLT ETRIG timestamps found! Check data quality!"; - } - else if (_ptb_hlt_etrig.size() == 1){ - _hlt_etrig = _ptb_hlt_etrig[0]; - _hlt_etrig_ts = _ptb_hlt_etrig_ts[0]; - } - else if(_ptb_hlt_etrig.size() > 1){ - uint64_t min_diff = std::numeric_limits::max(); - for(size_t i = 0; i < _ptb_hlt_etrig.size(); i++){ - uint64_t diff = _raw_ts > _ptb_hlt_etrig_ts[i] ? _raw_ts - _ptb_hlt_etrig_ts[i] : _ptb_hlt_etrig_ts[i] - _raw_ts; - if(diff < min_diff) - { - min_diff = diff; - _hlt_etrig = _ptb_hlt_etrig[i]; - _hlt_etrig_ts = _ptb_hlt_etrig_ts[i]; - } + else + { + uint64_t min_diff = std::numeric_limits::max(); + for(size_t i = 0; i < _ptb_hlt_etrig.size(); i++) + { + uint64_t diff = _raw_ts > _ptb_hlt_etrig_ts[i] ? _raw_ts - _ptb_hlt_etrig_ts[i] : _ptb_hlt_etrig_ts[i] - _raw_ts; + + if(diff < min_diff) + { + min_diff = diff; + _hlt_etrig = _ptb_hlt_etrig[i]; + _hlt_etrig_ts = _ptb_hlt_etrig_ts[i]; + } + } } - } +} - //Check which Gate/CRT T1 HLT to use based on ETRIG HLT - //Order to check: Beam -> Offbeam -> Xmuon - for (size_t i = 0; i < fBeamEtrigHlt.size(); i++){ - if (_hlt_etrig == fBeamEtrigHlt[i]){ - _hlt_gate = fBeamGateHlt; - _hlt_crtt1 = fBeamCrtT1Hlt; - _isBeam = true; - break; - } - } - if (!_isBeam){ - for (size_t i = 0; i < fOffbeamEtrigHlt.size(); i++){ - if (_hlt_etrig == fOffbeamEtrigHlt[i]){ - _hlt_gate = fOffbeamGateHlt; - _hlt_crtt1 = fOffbeamCrtT1Hlt; - _isOffbeam = true; - break; - } - } - } - //TODO: Grab Xmuon HLTs - - //Get Gate and CRT T1 HLT timestamps - //TODO: What if there is no Gate or CRT T1 HLT? - for (size_t i = 0; i < _ptb_hlt_unmask_timestamp.size(); i++){ - if(_ptb_hlt_trunmask[i] == _hlt_gate){ - _hlt_gate_ts = _ptb_hlt_unmask_timestamp[i]; +uint64_t sbnd::timing::FrameShift::DecideGlobalEtrigTimestamp() +{ + //Decide which global frame to use as reference + // Prioritise TDC ETRIG then PTB ETRIG then the raw header. + uint64_t global_etrig_ts = kInvalidTimestamp; + + if(_tdc_etrig_ts != kInvalidTimestamp) + global_etrig_ts = _tdc_etrig_ts; + else if(_hlt_etrig_ts != kInvalidTimestamp) + global_etrig_ts = _hlt_etrig_ts; + else + global_etrig_ts = _raw_ts; + + if(fDebugFrame) + { + std::cout << "----------------------------------------------------" << std::endl; + if(_tdc_etrig_ts != kInvalidTimestamp) + std::cout << "Using TDC ETRIG as Global Frame Reference" << std::endl; + else if(_hlt_etrig_ts != kInvalidTimestamp) + std::cout << "Using PTB HLT ETRIG as Global Frame Reference" << std::endl; + else + std::cout << "Using DAQ Header Timestamp as Global Frame Reference" << std::endl; + + std::cout << "Global Frame Timestamp: " << PrintFormatTimestamp(global_etrig_ts) << std::endl; } - if(_ptb_hlt_trunmask[i] == _hlt_crtt1){ - _hlt_crtt1_ts = _ptb_hlt_unmask_timestamp[i]; + + return global_etrig_ts; +} + +void sbnd::timing::FrameShift::DecideRelevantTDCTimestamps(const uint64_t &global_etrig_ts) +{ + // ch0: CRT T1 + if(_tdc_ch0.size() != 0) + _tdc_crtt1_ts = FindClosest(_tdc_ch0, global_etrig_ts); + + // ch1: BES + if(_tdc_ch1.size() != 0) + _tdc_bes_ts = FindClosest(_tdc_ch1, global_etrig_ts); + + // ch2: RWM + if(_tdc_ch2.size() != 0) + _tdc_rwm_ts = FindClosest(_tdc_ch2, global_etrig_ts); + + if(fDebugTdc) + { + std::cout << "----------------------------------------------------" << std::endl; + std::cout << "TDC Channel 0 (CRTT1) Timestamp: " << PrintFormatTimestamp(_tdc_crtt1_ts) << std::endl; + std::cout << "TDC Channel 1 (BES) Timestamp: " << PrintFormatTimestamp(_tdc_bes_ts) << std::endl; + std::cout << "TDC Channel 2 (RWM) Timestamp: "<< PrintFormatTimestamp(_tdc_rwm_ts) << std::endl; + std::cout << "TDC Channel 4 (ETRIG) Timestamp: "<< PrintFormatTimestamp(_tdc_etrig_ts) << std::endl; + std::cout << "----------------------------------------------------" << std::endl; + //for (size_t i = 0; i < _tdc_ch3.size(); i++) + // { + // std::cout << "TDC Channel 3 (FTRIG) Timestamp " << i << ": " << PrintFormatTimestamp(_tdc_ch3[i]) << std::endl; + // } + //std::cout << "----------------------------------------------------" << std::endl; } - } - - if (fDebugPtb){ - - std::cout << "----------------------------------------------------" << std::endl; - if (_isBeam) std::cout << "This is Beam Stream!" << std::endl; - if (_isOffbeam) std::cout << "This is Offbeam Stream!" << std::endl; - std::cout << "HLT ETRIG = " << _hlt_etrig - << ", Timestamp: " - << " (s) = " << _hlt_etrig_ts/uint64_t(1e9) - << ", (ns) = " << _hlt_etrig_ts%uint64_t(1e9) - << std::endl; - std::cout << "HLT Gate = " << _hlt_gate - << ", Timestamp: " - << " (s) = " << _hlt_gate_ts/uint64_t(1e9) - << ", (ns) = " << _hlt_gate_ts%uint64_t(1e9) - << std::endl; - std::cout << "HLT CRT T1 = " << _hlt_crtt1 - << ", Timestamp: " - << " (s) = " << _hlt_crtt1_ts/uint64_t(1e9) - << ", (ns) = " << _hlt_crtt1_ts%uint64_t(1e9) - << std::endl; - std::cout << "----------------------------------------------------" << std::endl; - } - - //--------------------------Get the Global Frame--------------------------// - - if (_global_timing_type == 0){ - _global_frame = _tdc_etrig_ts; - }else if(_global_timing_type ==1){ - _global_frame = _hlt_etrig_ts; - } - //TODO: Add for global timing type 2 (nothing) - - uint64_t diff = _raw_ts > _global_frame ? _raw_ts - _global_frame : _global_frame - _raw_ts; - if (diff > fMaxAllowedRefTimeDiff){ - throw cet::exception("FrameShift") << "Global timestamp is too far from the raw timestamp! " - << "Raw TS = (s) " << _raw_ts/uint64_t(1e9) << " (ns) " << _raw_ts%uint64_t(1e9) << std::endl - << ", Gobal Frame TS = (s) " << _global_frame/uint64_t(1e9) << " (ns) " << _global_frame%uint64_t(1e9) - << ", Diff = " << diff - << ". Check data quality!"; - } - if (fDebugFrame){ - std::cout << "----------------------------------------------------" << std::endl; - if(_global_timing_type == 0){ - std::cout << "Global Frame Type = SPECTDC (0)" << std::endl; - }else if(_global_timing_type == 1){ - std::cout << "Global Frame Type = PTB HLT (1)" << std::endl; - }else{ - std::cout << "Global Frame Type = Nothing (2)" << std::endl; +} + +void sbnd::timing::FrameShift::DecideRelevantPTBTimestamps(const uint64_t &global_etrig_ts) +{ + //Check which Gate/CRT T1 HLT to use based on ETRIG HLT + //Order to check: Beam -> Offbeam -> Xmuon + + for(const int &beam_etrig_hlt : fBeamEtrigHlt) + { + if(_hlt_etrig == beam_etrig_hlt) + { + _hlt_gate = fBeamGateHlt; + _hlt_crtt1 = fBeamCrtT1Hlt; + _isBeam = true; + break; + } } - std::cout << "Global Frame Timestamp: " - << " (s) = " << _global_frame/uint64_t(1e9) - << ", (ns) = " << _global_frame%uint64_t(1e9) - << std::endl; - std::cout << "----------------------------------------------------" << std::endl; - } - - //-----------------------Compute Frame Shift-----------------------// - - //for beam stream - if (_isBeam){ - if(_tdc_crtt1_ts != 0) _frame_tdc_crtt1 = SubtractUTCTimestmap(_global_frame, _tdc_crtt1_ts); - if(_tdc_bes_ts != 0) _frame_tdc_bes = SubtractUTCTimestmap(_global_frame, _tdc_bes_ts); - if(_tdc_rwm_ts != 0) _frame_tdc_rwm = SubtractUTCTimestmap(_global_frame, _tdc_rwm_ts); - if(_hlt_crtt1_ts != 0) _frame_hlt_crtt1 = SubtractUTCTimestmap(_global_frame, _hlt_crtt1_ts); - if(_hlt_gate_ts != 0) _frame_hlt_gate = SubtractUTCTimestmap(_global_frame, _hlt_gate_ts); - } - //for offbeam stream - else if (_isOffbeam){ - if(_tdc_crtt1_ts != 0) _frame_tdc_crtt1 = SubtractUTCTimestmap(_global_frame, _tdc_crtt1_ts); - if(_hlt_crtt1_ts != 0) _frame_hlt_crtt1 = SubtractUTCTimestmap(_global_frame, _hlt_crtt1_ts); - if(_hlt_gate_ts != 0) _frame_hlt_gate = SubtractUTCTimestmap(_global_frame, _hlt_gate_ts); - } - //TODO: for cosmic stream - - if (fDebugFrame){ - std::cout << "--------------------------------------" << std::endl; - std::cout << "Frame TDC CRT T1 = " << std::setprecision(9) << _frame_tdc_crtt1 << " (ns)" << std::endl; - std::cout << "Frame TDC BES = " << std::setprecision(9) << _frame_tdc_bes << " (ns)" << std::endl; - std::cout << "Frame TDC RWM = " << std::setprecision(9) << _frame_tdc_rwm << " (ns)" << std::endl; - std::cout << "--------------------------------------" << std::endl; - std::cout << "Frame HLT CRT T1 = " << std::setprecision(9) << _frame_hlt_crtt1 << " (ns)" << std::endl; - std::cout << "Frame HLT Gate = " << std::setprecision(9) << _frame_hlt_gate << " (ns)" << std::endl; - std::cout << "--------------------------------------" << std::endl; - } - - //-----------------------Pick which one to use in CAF-----------------------// - // The follow picks which frame to apply at the CAF stage and store it as frameApplyAtCaf, based on the stream - // This is so that no decision making is applied CAF. - // - // 1. Beam Stream: recontruct the beam spill + porch relative to RWM record in TDC since it has better resolution. - // Then, apply a constant to shift the RWM fram - // i.e. frame_tdc_rwm + shift rwm to beam gate HLT 26 - // - // 2. Offbeam Stream: reconstruct the porch relative to the beam gate opening frame - // i.e. frame_hlt_gate equiv - // - // 3. Xmuon Stream: do nothing for now - // - // TODO: Align Data Beam and Offbeam with MC with frame_data2mc - // + MC: t = 0 = first proton in spill - // + Data: t = 0 = abitrary. All subsystem electronics time is reference to the last PPS - - if(_isBeam){ - if (_frame_tdc_rwm != 0){ - _frame_tdc_rwm += fShiftRWM2Gate; //align RWM signal to gate opening frame - _frame_apply_at_caf = _frame_tdc_rwm; // +frame_data2mc - }else if(_frame_hlt_gate != 0){ - _frame_apply_at_caf = _frame_hlt_gate; //+frame_data2mc - }else{ - _frame_apply_at_caf = 0; + + if(!_isBeam) + { + for(const int &offbeam_etrig_hlt : fOffbeamEtrigHlt) + { + if(_hlt_etrig == offbeam_etrig_hlt) + { + _hlt_gate = fOffbeamGateHlt; + _hlt_crtt1 = fOffbeamCrtT1Hlt; + _isOffbeam = true; + break; + } + } } - } - else if(_isOffbeam){ - if(_frame_hlt_gate != 0){ - _frame_apply_at_caf = _frame_hlt_gate; //+frame_data2mc - }else{ - _frame_apply_at_caf = 0; + + if(!_isBeam & !_isOffbeam) + { + for(const int &xmuon_etrig_hlt : fXmuonEtrigHlt) + { + if(_hlt_etrig == xmuon_etrig_hlt) + { + _isXmuon = true; + break; + } + } } - } - else if(_isXmuon){ - _frame_apply_at_caf = 0; - } - //Put product in event - std::unique_ptr< sbnd::timing::FrameShiftInfo > newFrameShiftInfo(new sbnd::timing::FrameShiftInfo(_global_timing_type, _frame_tdc_crtt1, _frame_tdc_bes, _frame_tdc_rwm, _frame_hlt_crtt1, _frame_hlt_gate, _frame_apply_at_caf)); - std::unique_ptr< sbnd::timing::TimingInfo > newTimingInfo(new sbnd::timing::TimingInfo(_raw_ts, _tdc_crtt1_ts, _tdc_bes_ts, _tdc_rwm_ts, _tdc_etrig_ts, _hlt_crtt1_ts, _hlt_etrig_ts, _hlt_gate_ts)); + if( !_isBeam & !_isOffbeam & !_isXmuon) + throw cet::exception("FrameShift") << "ETRIG HLT " << _hlt_etrig << " does not match any known Beam/Offbeam/Xmuon ETRIG HLT! Check data quality!"; - e.put(std::move(newTimingInfo)); - e.put(std::move(newFrameShiftInfo)); + //Get Gate and CRT T1 HLT timestamps + for(size_t i = 0; i < _ptb_hlt_unmask_timestamp.size(); i++) + { + if(_ptb_hlt_trunmask[i] == _hlt_gate) + _hlt_gate_ts = _ptb_hlt_unmask_timestamp[i]; + if(_ptb_hlt_trunmask[i] == _hlt_crtt1) + _hlt_crtt1_ts = _ptb_hlt_unmask_timestamp[i]; + } - //Fill the tree - fTree->Fill(); + if(fDebugPtb) + { + std::cout << "----------------------------------------------------" << std::endl; + if(_isBeam) std::cout << "This is Beam Stream!" << std::endl; + if(_isOffbeam) std::cout << "This is Offbeam Stream!" << std::endl; + if(_isXmuon) std::cout << "This is Crossing Muon Stream!" << std::endl; + std::cout << "HLT ETRIG = " << _hlt_etrig << ", Timestamp: " << PrintFormatTimestamp(_hlt_etrig_ts) << std::endl; + std::cout << "HLT Gate = " << _hlt_gate << ", Timestamp: " << PrintFormatTimestamp(_hlt_gate_ts) << std::endl; + std::cout << "HLT CRT T1 = " << _hlt_crtt1 << ", Timestamp: " << PrintFormatTimestamp(_hlt_crtt1_ts) << std::endl; + std::cout << "----------------------------------------------------" << std::endl; + } } -void sbnd::FrameShift::beginJob() +void sbnd::timing::FrameShift::beginJob() { - // Implementation of optional member function here. - //Event Tree - fTree = tfs->make("events", ""); + if(fMakeTree) + { + fTree = tfs->make("events", ""); - fTree->Branch("run", &_run); - fTree->Branch("subrun", &_subrun); - fTree->Branch("event", &_event); + fTree->Branch("run", &_run); + fTree->Branch("subrun", &_subrun); + fTree->Branch("event", &_event); - //TDC - fTree->Branch("tdc_ch0", &_tdc_ch0); - fTree->Branch("tdc_ch1", &_tdc_ch1); - fTree->Branch("tdc_ch2", &_tdc_ch2); - fTree->Branch("tdc_ch4", &_tdc_ch4); - - fTree->Branch("tdc_crtt1_ts", &_tdc_crtt1_ts); - fTree->Branch("tdc_bes_ts", &_tdc_bes_ts); - fTree->Branch("tdc_rwm_ts", &_tdc_rwm_ts); - fTree->Branch("tdc_etrig_ts", &_tdc_etrig_ts); - - //PTB - //fTree->Branch("ptb_hlt_trigger", &_ptb_hlt_trigger); - //fTree->Branch("ptb_hlt_timestamp", &_ptb_hlt_timestamp); - fTree->Branch("ptb_hlt_unmask_timestamp", &_ptb_hlt_unmask_timestamp); - fTree->Branch("ptb_hlt_trunmask", &_ptb_hlt_trunmask); - - fTree->Branch("hlt_etrig", &_hlt_etrig); - fTree->Branch("hlt_etrig_ts", &_hlt_etrig_ts); - fTree->Branch("hlt_gate", &_hlt_gate); - fTree->Branch("hlt_gate_ts", &_hlt_gate_ts); - fTree->Branch("hlt_crtt1", &_hlt_crtt1); - fTree->Branch("hlt_crtt1_ts", &_hlt_crtt1_ts); - - //Timing Reference - fTree->Branch("crt_timing_ch", &_crt_timing_ch); - fTree->Branch("crt_timing_type", &_crt_timing_type); - fTree->Branch("pmt_timing_ch", &_pmt_timing_ch); - fTree->Branch("pmt_timing_type", &_pmt_timing_type); - - //Frame Shift - fTree->Branch("global_timing_type", &_global_timing_type); - fTree->Branch("global_frame", &_global_frame); - fTree->Branch("frame_tdc_crtt1", &_frame_tdc_crtt1); - fTree->Branch("frame_tdc_bes", &_frame_tdc_bes); - fTree->Branch("frame_tdc_rwm", &_frame_tdc_rwm); - fTree->Branch("frame_hlt_crtt1", &_frame_hlt_crtt1); - fTree->Branch("frame_hlt_gate", &_frame_hlt_gate); - fTree->Branch("frame_data2mc", &_frame_data2mc); - fTree->Branch("frame_apply_at_caf", &_frame_apply_at_caf); + //TDC + fTree->Branch("tdc_ch0", &_tdc_ch0); + fTree->Branch("tdc_ch1", &_tdc_ch1); + fTree->Branch("tdc_ch2", &_tdc_ch2); + fTree->Branch("tdc_ch4", &_tdc_ch4); + + fTree->Branch("tdc_crtt1_ts", &_tdc_crtt1_ts); + fTree->Branch("tdc_bes_ts", &_tdc_bes_ts); + fTree->Branch("tdc_rwm_ts", &_tdc_rwm_ts); + fTree->Branch("tdc_etrig_ts", &_tdc_etrig_ts); + + //PTB + fTree->Branch("ptb_hlt_unmask_timestamp", &_ptb_hlt_unmask_timestamp); + fTree->Branch("ptb_hlt_trunmask", &_ptb_hlt_trunmask); + + fTree->Branch("hlt_etrig", &_hlt_etrig); + fTree->Branch("hlt_etrig_ts", &_hlt_etrig_ts); + fTree->Branch("hlt_gate", &_hlt_gate); + fTree->Branch("hlt_gate_ts", &_hlt_gate_ts); + fTree->Branch("hlt_crtt1", &_hlt_crtt1); + fTree->Branch("hlt_crtt1_ts", &_hlt_crtt1_ts); + + //Frame Shift + fTree->Branch("frame_crtt1", &_frame_crtt1); + fTree->Branch("timing_type_crtt1", &_timing_type_crtt1); + fTree->Branch("timing_channel_crtt1", &_timing_channel_crtt1); + + fTree->Branch("frame_gate", &_frame_gate); + fTree->Branch("timing_type_gate", &_timing_type_gate); + fTree->Branch("timing_channel_gate", &_timing_channel_gate); + + fTree->Branch("frame_etrig", &_frame_etrig); + fTree->Branch("timing_type_etrig", &_timing_type_etrig); + fTree->Branch("timing_channel_etrig", &_timing_channel_etrig); + + fTree->Branch("frame_default", &_frame_default); + fTree->Branch("timing_type_default", &_timing_type_default); + fTree->Branch("timing_channel_default", &_timing_channel_default); + } } -void sbnd::FrameShift::ResetEventVars() +void sbnd::timing::FrameShift::ResetEventVars() { _run = -1; _subrun = -1; _event = -1; - _raw_ts = 0; + _raw_ts = kInvalidTimestamp; _tdc_ch0.clear(); _tdc_ch1.clear(); _tdc_ch2.clear(); + //_tdc_ch3.clear(); _tdc_ch4.clear(); - _tdc_crtt1_ts = 0; - _tdc_bes_ts = 0; - _tdc_rwm_ts = 0; - _tdc_etrig_ts = 0; + _tdc_crtt1_ts = kInvalidTimestamp; + _tdc_bes_ts = kInvalidTimestamp; + _tdc_rwm_ts = kInvalidTimestamp; + _tdc_etrig_ts = kInvalidTimestamp; _ptb_hlt_trigger.clear(); _ptb_hlt_timestamp.clear(); @@ -773,57 +741,29 @@ void sbnd::FrameShift::ResetEventVars() _isBeam = false; _isOffbeam = false; _isXmuon = false; - _hlt_etrig = -1; - _hlt_etrig_ts = 0; - _hlt_gate = -1; - _hlt_gate_ts = 0; - _hlt_crtt1= -1; - _hlt_crtt1_ts = 0; - - _crt_timing_ch = 99; - _crt_timing_type = 99; - _pmt_timing_ch = 99; - _pmt_timing_type = 99; - - _global_timing_type = 99; - _global_frame = 0; - - _frame_tdc_crtt1 = 0; - _frame_tdc_bes = 0; - _frame_tdc_rwm = 0; - _frame_hlt_crtt1 = 0; - _frame_hlt_gate = 0; - _frame_apply_at_caf = 0; -} -double sbnd::FrameShift::SubtractUTCTimestmap(const uint64_t& ts1, const uint64_t& ts2) -{ - //Subtract two timestamps in UTC format - //ts1 and ts2 are in nanoseconds in uint64_t format - //Return the difference in nanoseconds as double - - double ts1_s = ts1 / uint64_t(1e9); - double ts1_ns = ts1 % uint64_t(1e9); - double ts2_s = ts2 / uint64_t(1e9); - double ts2_ns = ts2 % uint64_t(1e9); - - //std::cout << std::setprecision(15) << "ts1_s = " << ts1_s << std::setprecision(9)<< ", ts1_ns = " << ts1_ns << std::endl; - //std::cout << std::setprecision(15) << "ts2_s = " << ts2_s << std::setprecision(9)<< ", ts2_ns = " << ts2_ns << std::endl; + _hlt_etrig = std::numeric_limits::max(); + _hlt_etrig_ts = kInvalidTimestamp; + _hlt_gate = std::numeric_limits::max(); + _hlt_gate_ts = kInvalidTimestamp; + _hlt_crtt1= std::numeric_limits::max(); + _hlt_crtt1_ts = kInvalidTimestamp; - double diff_s = 0; - double diff_ns = 0; + _frame_crtt1 = kInvalidTimestamp; + _timing_type_crtt1 = kInvalidType; + _timing_channel_crtt1 = kInvalidChannel; - //If the same PPS, just subtract the nanoseconds - if(ts1_s == ts2_s){ - diff_ns = ts1_ns - ts2_ns; - } - //If ts1 is later than ts2, then subtract the seconds and add the nanoseconds - else{ - diff_s = ts1_s - ts2_s; - diff_ns = diff_s * (double)1e9 + ts1_ns - ts2_ns; - } - - return diff_ns; + _frame_gate = kInvalidTimestamp; + _timing_type_gate = kInvalidType; + _timing_channel_gate = kInvalidChannel; + + _frame_etrig = kInvalidTimestamp; + _timing_type_etrig = kInvalidType; + _timing_channel_etrig = kInvalidChannel; + + _frame_default = kInvalidTimestamp; + _timing_type_default = kInvalidType; + _timing_channel_default = kInvalidChannel; } -DEFINE_ART_MODULE(sbnd::FrameShift) +DEFINE_ART_MODULE(sbnd::timing::FrameShift) diff --git a/sbndcode/Timing/FrameShift/frameshift_sbnd_data.fcl b/sbndcode/Timing/FrameShift/frameshift_sbnd_data.fcl index 739fbc792..6abc07e1f 100644 --- a/sbndcode/Timing/FrameShift/frameshift_sbnd_data.fcl +++ b/sbndcode/Timing/FrameShift/frameshift_sbnd_data.fcl @@ -1,40 +1,40 @@ BEGIN_PROLOG + frameshift_data: { - module_type: "FrameShift" - - DAQHeaderInstanceLabel: "RawEventHeader" - DAQHeaderModuleLabel: "daq" - RawTSCorrection: 367000 //367 us - MaxAllowedRefTimeDiff: 3000000 //3 ms + module_type: "FrameShift" - TimingRefPmtLabel: "pmtdecoder" - TimingRefCrtLabel: "crtstrips" - - TdcDecodeLabel: "tdcdecoder" - PtbDecodeLabel: "ptbdecoder" + DAQHeaderInstanceLabel: "RawEventHeader" + DAQHeaderModuleLabel: "daq" + RawTSCorrection: 367000 #367 us + + TdcDecodeLabel: "tdcdecoder" + PtbDecodeLabel: "ptbdecoder" - PtbEtrigHlts: [1, 2, 3, 4, 5, 14, 15] - - BeamEtrigHlt: [1,2] - OffbeamEtrigHlt: [3,4] - XmuonEtrigHlt: [5, 14, 15] + PtbEtrigHlts: [1, 2, 3, 4, 5, 14, 15, 16, 17] + + BeamEtrigHlt: [1,2,16] + OffbeamEtrigHlt: [3,4,17] + XmuonEtrigHlt: [5, 14, 15] - BeamCrtT1Hlt: 20 - OffbeamCrtT1Hlt: 21 + BeamCrtT1Hlt: 20 + OffbeamCrtT1Hlt: 21 - BeamGateHlt: 26 - OffbeamGateHlt: 27 + BeamGateHlt: 26 + OffbeamGateHlt: 27 - DebugDAQHeader: false - DebugTimingRef: false - DebugPtb: false - DebugTdc: false - DebugFrame: false - - ShiftData2MC: 0 - ShiftRWM2Gate: 1738 -} - + DebugDAQHeader: false + DebugPtb: false + DebugTdc: false + DebugFrame: false + ShiftData2MC: 0 + ShiftTdcRwm2PtbGate: 1738 #ns + #TDC RWM arrives 1738 ns later than HLT Gate on average, so shift TDC RWM earlier + ShiftTdcEtrig2PtbEtrig: 257 #ns + #TDC ETRIG arrives 257 ns later than HLT ETRIG on average, so shift TDC ETRIG earlier + ShiftTdcCrtt12PtbCrtt1: 295 #ns + #TDC CRTT1 arrives 295 ns later than HLT CRTT1 on average, so shift TDC CRTT1 earlier +} + END_PROLOG diff --git a/sbndcode/Timing/TimingUtils.cc b/sbndcode/Timing/TimingUtils.cc new file mode 100644 index 000000000..06ec201e1 --- /dev/null +++ b/sbndcode/Timing/TimingUtils.cc @@ -0,0 +1,27 @@ +#include "sbndcode/Timing/TimingUtils.h" + +namespace sbnd { + double TimingUtils::SubtractUTCTimestmap(const uint64_t& ts1, const uint64_t& ts2) + { + + double ts1_s = ts1 / uint64_t(1e9); + double ts1_ns = ts1 % uint64_t(1e9); + double ts2_s = ts2 / uint64_t(1e9); + double ts2_ns = ts2 % uint64_t(1e9); + + double diff_s = 0; + double diff_ns = 0; + + //If the same PPS, just subtract the nanoseconds + if(ts1_s == ts2_s){ + diff_ns = ts1_ns - ts2_ns; + } + //If ts1 is later than ts2, then subtract the seconds and add the nanoseconds + else{ + diff_s = ts1_s - ts2_s; + diff_ns = diff_s * (double)1e9 + ts1_ns - ts2_ns; + } + + return diff_ns; + } +} diff --git a/sbndcode/Timing/TimingUtils.h b/sbndcode/Timing/TimingUtils.h new file mode 100644 index 000000000..6b84cbf2b --- /dev/null +++ b/sbndcode/Timing/TimingUtils.h @@ -0,0 +1,24 @@ +#ifndef TIMINGUTILS_H_SEEN +#define TIMINGUTILS_H_SEEN + +/////////////////////////////////////////////// +// TimingUtils.h +// +// Common functions for Timing reconstruction +/////////////////////////////////////////////// + +#include "canvas/Persistency/Common/Ptr.h" + +namespace sbnd { + + namespace TimingUtils{ + + //Subtract two timestamps in UTC format + //ts1 and ts2 are in nanoseconds in uint64_t format + //Return the difference in nanoseconds as double + double SubtractUTCTimestmap(const uint64_t& ts1, const uint64_t& ts2); + + } +} + +#endif \ No newline at end of file diff --git a/sbndcode/Timing/WaveformAlignment/WaveformAlignment_module.cc b/sbndcode/Timing/WaveformAlignment/WaveformAlignment_module.cc index d041b580f..636e9004c 100644 --- a/sbndcode/Timing/WaveformAlignment/WaveformAlignment_module.cc +++ b/sbndcode/Timing/WaveformAlignment/WaveformAlignment_module.cc @@ -39,6 +39,7 @@ #include "lardataobj/RawData/OpDetWaveform.h" #include "sbnobj/SBND/Timing/DAQTimestamp.hh" +#include "sbnobj/SBND/Timing/FrameShiftInfo.hh" #include "sbndcode/Timing/SBNDRawTimingObj.h" #include "sbndcode/Decoders/PTB/sbndptb.h" #include "sbndcode/Calibration/PDSDatabaseInterface/PMTCalibrationDatabase.h" @@ -143,7 +144,6 @@ class sbnd::WaveformAlignment : public art::EDProducer { // TDC stuff std::vector _tdc_ch3; std::vector _tdc_ch3_ps; - std::vector _tdc_ch4; //PTB stuff std::vector _ptb_hlt_trigger; @@ -157,7 +157,7 @@ class sbnd::WaveformAlignment : public art::EDProducer { //---FHICL CONFIG PARAMETERS // Product label - art::InputTag fTimingRefLabel; + art::InputTag fFrameShiftLabel; art::InputTag fTdcDecodeLabel; art::InputTag fPtbDecodeLabel; art::InputTag fPmtDecodeLabel; @@ -170,6 +170,9 @@ class sbnd::WaveformAlignment : public art::EDProducer { //PTB std::vector fPtbAllowedHlts; + //FrameShift + uint64_t _frameDefault; + // DAQ double fWfLength; int fnChperBoard; @@ -233,7 +236,7 @@ sbnd::WaveformAlignment::WaveformAlignment(fhicl::ParameterSet const& p) : EDProducer{p} // // More initializers here. { - fTimingRefLabel = p.get("TimingRefLabel", "pmtdecoder"); + fFrameShiftLabel = p.get("FrameShiftLabel", "frameshift"); fTdcDecodeLabel = p.get("TdcDecodeLabel", "tdcdecoder"); fPtbDecodeLabel = p.get("PtbDecodeLabel", "ptbdecoder"); @@ -368,27 +371,18 @@ void sbnd::WaveformAlignment::produce(art::Event& e) if (fDebugTimeRef | fDebugFtrig | fDebugPmt | fDebugTiming) std::cout <<"#----------RUN " << _run << " SUBRUN " << _subrun << " EVENT " << _event <<"----------#\n"; - //------------------------PMT Timing--------------------------// - art::Handle timingRefHandle; - e.getByLabel(fTimingRefLabel, timingRefHandle); + //------------------------Frame Shift--------------------------// + art::Handle frameShiftHandle; + e.getByLabel(fFrameShiftLabel, frameShiftHandle); - if (!timingRefHandle.isValid()){ - throw cet::exception("WaveformAlignment") << "No raw::TimingReferenceInfo found w/ tag " << fTimingRefLabel << ". Check data quality!"; + if (!frameShiftHandle.isValid()){ + throw cet::exception("WaveformAlignment") << "No sbnd::timing::FrameShiftInfo found w/ tag " << fFrameShiftLabel << ". Check data quality!"; } else{ - raw::TimingReferenceInfo const& pmt_timing(*timingRefHandle); - - _pmt_timing_type = pmt_timing.timingType; - _pmt_timing_ch = pmt_timing.timingChannel; - - if (fDebugTimeRef){ - std::cout << "Timing Reference For Decoding PMT" << std::endl; - std::cout << " Type = " << _pmt_timing_type << " (SPECTDC = 0; PTB HLT = 1; CAEN-only = 3)." << std::endl; - std::cout << " Channel = " << _pmt_timing_ch << " (TDC ETRIG = 4; PTB BNB Beam+Light = 2)." << std::endl; - } + _frameDefault = frameShiftHandle->FrameDefault(); + _pmt_timing_type = frameShiftHandle->TimingTypeDefault(); } - timingRefHandle.removeProduct(); - + //---------------------------TDC-----------------------------// art::Handle> tdcHandle; e.getByLabel(fTdcDecodeLabel, tdcHandle); @@ -408,12 +402,9 @@ void sbnd::WaveformAlignment::produce(art::Event& e) const uint32_t ch = tdc.Channel(); const uint64_t ts = tdc.Timestamp(); - if ((ch == 3) || (ch ==4)){ - if (ch == 3) { - _tdc_ch3.push_back(ts); - _tdc_ch3_ps.push_back(tdc.TimestampPs()); - } - if (ch == 4) _tdc_ch4.push_back(ts); + if (ch == 3){ + _tdc_ch3.push_back(ts); + _tdc_ch3_ps.push_back(tdc.TimestampPs()); if (fDebugTimeRef) std::cout << "Ch " << name << " sec (s) = " << ts/uint64_t(1e9) @@ -560,21 +551,21 @@ void sbnd::WaveformAlignment::produce(art::Event& e) }else{ //TDC as reference frame if timing type is 0 if (_pmt_timing_type == 0){ - if (_tdc_ch4.size() == 0 || _tdc_ch3.size() == 0 || _tdc_ch3_ps.size() == 0){ + if ( _tdc_ch3.size() == 0 || _tdc_ch3_ps.size() == 0){ status = kMissingTDC; }else{ double tempShift = 0; if (_tdc_ch3_ps.front() == 0){ //old decode version does not have picosecond - tempShift = ComputeShift( _tdc_ch4.front() - , wf->TimeStamp() + tempShift = ComputeShift( _frameDefault //decodeFrame + , wf->TimeStamp() //decodeTs , tickFtrig , _tdc_ch3 , fTdc3CaenOffset[wf->ChannelNumber()] , false ); }else{ - tempShift = ComputeShiftPs( _tdc_ch4.front() - , wf->TimeStamp() + tempShift = ComputeShiftPs( _frameDefault //decodeFrame + , wf->TimeStamp() //decodeTs , tickFtrig , _tdc_ch3 , _tdc_ch3_ps @@ -597,20 +588,10 @@ void sbnd::WaveformAlignment::produce(art::Event& e) //OR timing type is 0 but cannot find the TDC frame //TDC status error will be replaced with PTB status error if((_pmt_timing_type == 1) | (status == kMissingTDC) | (status == kOutOfBoundTDC)){ - //Find the HLT of ETRIG - uint64_t tsFrame = std::numeric_limits::min(); - for (size_t i = 0; i < _ptb_hlt_unmask_timestamp.size(); i++){ - for (size_t j = 0; j < fPtbAllowedHlts.size(); j++){ - if (_ptb_hlt_trunmask[i] == fPtbAllowedHlts[j]){ - tsFrame = _ptb_hlt_unmask_timestamp[i]; - break; - } - } - } - if(tsFrame == std::numeric_limits::min() || _ptb_hlt_unmask_timestamp.size() == 0){ + if( _ptb_hlt_unmask_timestamp.size() == 0){ status = kMissingPTB; }else{ - double tempShift = ComputeShift( tsFrame + double tempShift = ComputeShift( _frameDefault , wf->TimeStamp() , tickFtrig , _ptb_hlt_unmask_timestamp @@ -877,7 +858,6 @@ void sbnd::WaveformAlignment::beginJob() fTree->Branch("event", &_event); fTree->Branch("tdc_ch3", &_tdc_ch3); fTree->Branch("tdc_ch3_ps", &_tdc_ch3_ps); - fTree->Branch("tdc_ch4", &_tdc_ch4); fTree->Branch("pmt_timing_type", &_pmt_timing_type); fTree->Branch("pmt_timing_ch", &_pmt_timing_ch); @@ -911,12 +891,13 @@ void sbnd::WaveformAlignment::ResetEventVars() { _run = -1; _subrun = -1; _event = -1; + _frameDefault = sbnd::timing::kInvalidTimestamp; + _pmt_timing_type = std::numeric_limits::max(); _pmt_timing_ch = std::numeric_limits::max(); _tdc_ch3.clear(); _tdc_ch3_ps.clear(); - _tdc_ch4.clear(); _ptb_hlt_trigger.clear(); _ptb_hlt_timestamp.clear(); @@ -1168,7 +1149,7 @@ bool sbnd::WaveformAlignment::CheckShift(const double shift, const int boardId) if (boardId == 8){ if ((shift > fTimingJitterUpperBound) | (shift < fTimingJitterLowerBound)){ if (fDebugFtrig) std::cout << " board id = " << boardId - << " has jittering " << shift + << " has jittering " << std::setprecision(9) << shift << " out of bound [" << fTimingJitterLowerBound << ", " << fTimingJitterUpperBound << "]" << std::endl; return false; @@ -1176,7 +1157,7 @@ bool sbnd::WaveformAlignment::CheckShift(const double shift, const int boardId) }else{ if ((shift > fPmtJitterUpperBound) | (shift < fPmtJitterLowerBound)){ if (fDebugFtrig) std::cout << " board id = " << boardId - << " has jittering " << shift + << " has jittering " << std::setprecision(9)<< shift << " out of bound [" << fPmtJitterLowerBound << ", " << fPmtJitterUpperBound << "]" << std::endl; return false; diff --git a/sbndcode/Timing/WaveformAlignment/wfalign_sbnd_data.fcl b/sbndcode/Timing/WaveformAlignment/wfalign_sbnd_data.fcl index 1cf6d4d96..87bfd1f5d 100644 --- a/sbndcode/Timing/WaveformAlignment/wfalign_sbnd_data.fcl +++ b/sbndcode/Timing/WaveformAlignment/wfalign_sbnd_data.fcl @@ -20,7 +20,7 @@ wfalign_data: #Configure input projects #Types of High Level Trigger for data stream - PtbAllowedHlts: [1, 2, 3, 4, 5, 14, 15] + PtbAllowedHlts: [1, 2, 3, 4, 5, 14, 15, 16, 17] #Time interval between TDC ch3 (FTRIG) and ch15 CAEN Tdc3CaenOffset: [101, 101, 101, 101, 101, 101, 101, 101, 101] #ns #Time interval between HLT issued in PTB and ch15 CAEN