EliminateBeamSpills_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief A new version of RemoveBeamSpills that takes a list of
3 /// NuMI event times as input, and does not depend on the
4 /// mysterious database that RemoveBeamSpills does.
5 ///
6 /// You must generate a list of times using SpillTime_module.cc. These
7 /// times are NOvA times, a monotonic timescale. See more details in
8 /// SpillTime. The file format is a list of ASCII numbers separated by
9 /// whitespace.
10 ///
11 /// You can provide a file that extends considerably beyond
12 /// the bounds of the data you want to filter. The only significant harm
13 /// in providing an overly large list is memory usage. Each year of NuMI
14 /// triggers uses about 170MB of RAM. There is also a, probably negligible,
15 /// cost of O(n log n) to sort the list in the first place and O(log n)
16 /// to search in the list for each event.
17 ///
18 /// \author M. Strait <mstrait@fnal.gov>
19 ////////////////////////////////////////////////////////////////////////
20 
23 #include "DAQDataFormats/RawEvent.h"
24 #include "DAQDataFormats/RawDataBlock.h"
25 #include "RawData/FlatDAQData.h"
26 #include "RawData/RawTrigger.h"
27 
28 #include <string>
29 #include <algorithm>
30 #include <fstream>
31 
32 const int TDC_PER_US = 64;
33 
35 public:
36 
37  explicit EliminateBeamSpills(const fhicl::ParameterSet & pset): EDFilter(),
38  bufferbefore(pset.get<double>("bufferbefore")),
39  bufferafter(pset.get<double>("bufferafter")),
40  spillfile(pset.get<std::string>("spillfile")) { };
41 
42  virtual ~EliminateBeamSpills() { };
43  bool filter(art::Event& evt);
44  void beginJob();
45 
46 private:
47 
48  // If the trigger time of an event is before the beginning spill by
49  // less than this amount, in seconds, it will be removed by the filter.
50  double bufferbefore;
51 
52  // If the trigger time of an event is after the *beginning* of a
53  // spill (not the end, because it's not trivial to determine the
54  // length of each spill) by less than this amount, in seconds, it
55  // will be removed by the filter.
56  double bufferafter;
57 
58  // Name of the file containing the list of spill times. This file should
59  // contain a whitespace-delimited list of NOvA times of the spill (64ths of
60  // microseconds since the NOvA epoch).
62 
63  std::vector<uint64_t> spilltimes;
64 };
65 
66 
68 {
69  std::ifstream in(spillfile.c_str());
70  if(!in.is_open()){
71  LOG_ERROR("EliminateBeamSpills") << "Could not open spill timestamp file "
72  << spillfile;
73  abort();
74  }
75 
76  uint64_t t;
77 
78  while(in >> t) spilltimes.push_back(t);
79  spilltimes.shrink_to_fit();
80  std::sort(spilltimes.begin(), spilltimes.end());
81 }
82 
83 static bool start_and_end(uint64_t & event_start_time,
84  uint64_t & event_end_time,
85  const art::Handle< std::vector<rawdata::FlatDAQData> > & flatdaq,
86  const art::Handle< std::vector<rawdata::RawTrigger> > & rawtrigger)
87 {
89  if(flatdaq->empty()) return false;
90 
91  raw.readData((*flatdaq)[0].getRawBufferPointer());
92  if(raw.getDataBlockNumber() == 0) return false;
93 
94  raw.setFloatingDataBlock(0);
95  daqdataformats::RawDataBlock& datablock = *raw.getFloatingDataBlock();
96 
97  event_start_time = 0xffffffffffffffff;
98  event_end_time = 0x0000000000000000;
99 
100  for(unsigned int di = 0; di < raw.getDataBlockNumber(); di++){
101  raw.setFloatingDataBlock(di);
102  datablock = (*raw.getFloatingDataBlock());
103 
104  if(datablock.getHeader()->getMarker() ==
106  !datablock.getHeader()->checkMarker()) continue;
107 
108  for(unsigned int mi = 0; mi < datablock.getNumMicroBlocks(); mi++){
109  datablock.setFloatingMicroBlock(mi);
110  daqdataformats::RawMicroBlock * ub = datablock.getFloatingMicroBlock();
111 
112  // The time is always in the second and third words of the
113  // microslice, which follows two words of microblock header, so
114  // just get it. Justin says you can also get it from getTime(),
115  // but this already works and I'm not touching it.
116  const uint32_t t_marker_low = ((uint32_t *)(ub->getBuffer()))[3];
117  const uint32_t t_marker_high = ((uint32_t *)(ub->getBuffer()))[4];
118 
119  uint64_t time_marker = t_marker_low;
120  time_marker |= (uint64_t)t_marker_high << 32;
121  if(time_marker < event_start_time) event_start_time = time_marker;
122  if(time_marker > event_end_time ) event_end_time = time_marker;
123  }
124  }
125 
126  // Assume that microblocks are always 50us. I hope that's true for all
127  // relevant data.
128  event_end_time += 50 * TDC_PER_US;
129  return true; // ok
130 }
131 
132 static void getrawtrigger(
133  art::Handle< std::vector<rawdata::RawTrigger> > & trg,
134  const art::Event & evt)
135 {
136  evt.getByLabel("minbias", trg);
137  if(trg.failedToGet()) evt.getByLabel("daq", trg);
138 }
139 
140 static std::string prettyprintus(const double us)
141 {
142  char buf[256];
143  if (fabs(us) < 1e3) snprintf(buf, 256, "%.3f us", us);
144  else if(fabs(us) < 1e6) snprintf(buf, 256, "%.6f ms", us/1e3);
145  else snprintf(buf, 256, "%.9f s", us/1e6);
146  return buf;
147 }
148 
149 // Returns true if the event is to be kept, false if it is too close to a spill
151 {
152  if(spilltimes.empty()){
153  LOG_INFO("EliminateBeamSpills")
154  << "Empty spill time list. Hopefully on purpose? Passing.";
155  return true;
156  }
157 
159  evt.getByLabel("daq", flatdaq);
160 
162  getrawtrigger(rawtrigger, evt);
163  if(rawtrigger->empty()){
164  LOG_INFO("EliminateBeamSpills") << "No rawtrigger, failing event";
165  return false;
166  }
167 
168  uint64_t start, end;
169 
170  if(!start_and_end(start, end, flatdaq, rawtrigger)){
171  LOG_WARNING("EliminateBeamSpills")
172  << "Could not get event start and end, failing event";
173  return false;
174  }
175 
176  // Take all spills to be 10us long. Could be wrong depending on what AD is
177  // doing, but shouldn't be off by more than a few us. From doc-11535, I
178  // think the spills are slightly *less* than 10us, more like 9.8us, so this
179  // should be safe. (The user is encouraged to add a substantial buffer
180  // anyway).
181  const double spill_len_us = 10.0;
182  const uint64_t spill_len_tdc = spill_len_us * TDC_PER_US;
183 
184  const uint64_t startbuf = start - (uint64_t)(bufferbefore * TDC_PER_US * 1e6);
185  const uint64_t endbuf = end + (uint64_t)(bufferafter * TDC_PER_US * 1e6);
186 
187  // Find the next spill after the end of the event and its buffer
188  auto nextspill = upper_bound(spilltimes.begin(), spilltimes.end(), endbuf);
189 
190  // If all spills are after the event and its buffer, then we're good. But
191  // this may also signal that the user didn't provide all relevant spills.
192  if(nextspill == spilltimes.begin()){
193  LOG_WARNING("EliminateBeamSpills")
194  << "Next spill is "
195  << prettyprintus(double((int64_t)(*nextspill)-(int64_t)endbuf)/TDC_PER_US)
196  << " after the event and its buffer, but found no previous spill.\n"
197  "Passing with reservations.";
198  return true;
199  }
200 
201  if(nextspill == spilltimes.end()){
202  const auto previousspill = spilltimes.end()-1;
203  const uint64_t prevspillend = *previousspill + spill_len_tdc;
204  const bool ret = prevspillend <= startbuf;
205 
206  LOG_WARNING("EliminateBeamSpills") << (ret?"Passing":"Failing")
207  << ", end of previous spill is "
208  << prettyprintus(double((int64_t)startbuf-(int64_t)prevspillend)/TDC_PER_US)
209  << " from window start\nFound no next spill. Passing with reservations.";
210 
211  return ret;
212  }
213 
214  const auto previousspill = nextspill-1;
215  const uint64_t prevspillend = *previousspill + spill_len_tdc;
216  const bool ret = prevspillend <= startbuf;
217 
218  LOG_INFO("EliminateBeamSpills") << (ret?"Passing":"Failing")
219  << ", end of previous spill is "
220  << prettyprintus(double((int64_t)startbuf-(int64_t)prevspillend)/TDC_PER_US)
221  << " from window start\nAnd beginning of next spill is "
222  << prettyprintus(double((int64_t)(*nextspill)-(int64_t)endbuf)/TDC_PER_US)
223  << " from window end\n";
224 
225  return ret;
226 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
static std::string prettyprintus(const double us)
EliminateBeamSpills(const fhicl::ParameterSet &pset)
DEFINE_ART_MODULE(TestTMapFile)
const int TDC_PER_US
A new version of RemoveBeamSpills that takes a list of NuMI event times as input, and does not depend...
bool filter(art::Event &evt)
#define LOG_WARNING(category)
static constexpr double us
static const double ub
Definition: Units.h:88
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ifstream in
Definition: comparison.C:7
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
static bool start_and_end(uint64_t &event_start_time, uint64_t &event_end_time, const art::Handle< std::vector< rawdata::FlatDAQData > > &flatdaq, const art::Handle< std::vector< rawdata::RawTrigger > > &rawtrigger)
std::vector< uint64_t > spilltimes
#define LOG_INFO(stream)
Definition: Messenger.h:144
#define LOG_ERROR(stream)
Definition: Messenger.h:129
static void getrawtrigger(art::Handle< std::vector< rawdata::RawTrigger > > &trg, const art::Event &evt)