RemoveBeamSpills_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Reject triggers that are too close to beam spills
3 /// \author Christopher Backhouse - bckhouse@caltech.edu
4 ////////////////////////////////////////////////////////////////////////
5 
6 #include <string>
7 #include <iostream>
8 #include <time.h>
9 
16 #include "fhiclcpp/ParameterSet.h"
17 
18 #include "Geometry/Geometry.h"
19 #include "NovaTimingUtilities/TimingUtilities.h"
20 #include "RawData/RawTrigger.h"
21 
22 #include "TH1.h"
23 
24 #include "IFBeam_service.h"
25 
26 namespace util
27 {
28  /// Reject triggers that are too close to beam spills
30  {
31  public:
32  explicit RemoveBeamSpills(const fhicl::ParameterSet& pset);
33  virtual ~RemoveBeamSpills();
34 
35  void beginJob();
36  bool filter(art::Event& evt);
37 
38  protected:
41 
43 
44  std::unique_ptr<BeamFolder> fFolder;
45 
46  TH1* fClosest;
47  };
48 
49  //..............................................................
51  : fTotTrigs(0), fRejTrigs(0)
52  {
53  fRawDataLabel = pset.get<std::string>("RawDataLabel");
54  fVetoWindow = pset.get<double>("VetoWindow");
55  fEpsilon = pset.get<double>("Epsilon");
56 
58  fFolder = ifbeam->getBeamFolder(pset.get<std::string>("Bundle"),
59  pset.get<std::string>("URL"),
60  pset.get<double>("DBTimeWindow"));
61 
62  // Look this far back in time for a spill, in seconds
63  fFolder->set_epsilon(fEpsilon);
64  }
65 
66  //......................................................................
68  {
69  if(fTotTrigs)
70  std::cout << "RemoveBeamSpills: of " << fTotTrigs
71  << " triggers, rejected "
72  << fRejTrigs << " ("
73  << .01*((10000*fRejTrigs)/fTotTrigs) << "%)" << std::endl;
74  }
75 
76  //......................................................................
78  {
80  fClosest = tfs->make<TH1F>("closest", ";#Deltat (sec)", 2000, -.01, +.01);
81  }
82 
83  //......................................................................
85  {
86  // Doesn't apply to MC
87  if(!evt.isRealData()) return true;
88 
89  ++fTotTrigs;
90 
91  // Get trigger information to extract trigger time
93  evt.getByLabel(fRawDataLabel, trigv);
94  assert(!trigv->empty());
95  const rawdata::RawTrigger& trig = (*trigv)[0];
96 
97  // Get trigger time in unix time
98  struct timespec unixTime;
99  novadaq::timeutils::convertNovaTimeToUnixTime(trig.fTriggerTimingMarker_TimeStart, unixTime);
100  const double eventTime = unixTime.tv_sec + 1e-9*unixTime.tv_nsec;
101 
102  // If we're in the far detector, apply a time-of-flight correction to the trigger time
103  double tof = 0;
105  if (geom->DetId() == novadaq::cnv::kFARDET) tof = 810e3/3e8; // 810km/c
106  const double eventTimeCorrected = eventTime - tof;
107 
108  // Query the database to get NuMI trigger time. This query gives the EA9SNC value, which is
109  // a value stored to the database once per second that represents the time in seconds since
110  // the last trigger. Here we fetch both "ea9snc" (time since last trigger) and "time" (time
111  // value was stored to the database), and so "time - ea9snc" gives us the absolute time of
112  // the most recent trigger once per second.
113  //
114  // Note that we round the trigger time up to the closest second by calling ceil when querying
115  // the database below. This is intentional and makes ifbeam queries more robust. Querying the
116  // database can be a little fragile sometimes, and we can simplify it by making two assumptions:
117  //
118  // 1) the NuMI spill occurs at a rate less than 1Hz
119  // 2) Values are written to the database once per second, on the second
120  //
121  // If we round up the trigger time to the next second, then we'll pull out the value from the
122  // database that comes in at the _end_ of the second in which the trigger occurred. We then
123  // use the ea9snc value to look backwards in time to the most recent trigger; if the value is
124  // <1s then there's a chance it will coincide with the trigger and we veto it. The value stored
125  // in the database at the _start_ of the second in which the trigger occurs is not going to be
126  // useful to us, since the ea9snc value can only refer backwards in time to before our trigger
127  // occurred.
128 
129  double ea9snc, time;
130  try{
131  fFolder->GetNamedData(ceil(eventTimeCorrected), "G:EA9SNC@", &ea9snc, &time);
132  }
133  catch(std::exception& e) {
134  std::string err(e.what());
135  if (err.find("variable not found")) {
136  const time_t t = eventTimeCorrected;
137  mf::LogError("RemoveBeamSpills") << "Couldn't find anything in ifbeam db! Timestamp: "
138  << std::ctime(&t);
139  } else {
140  mf::LogError("RemoveBeamSpills") << "Exception encountered: " << e.what();
141  }
142  abort();
143  }
144  double deltaT = fabs(eventTimeCorrected - (time - ea9snc));
145  bool reject = deltaT < fVetoWindow;
146  if (reject) {
147  ++fRejTrigs;
148  mf::LogInfo("RemoveBeamSpills") << "Vetoed event with trigger time "
149  << std::fixed << eventTimeCorrected << " due to spill at time "
150  << time-ea9snc << " (deltaT = " << deltaT << ")" << std::defaultfloat;
151  }
152  return !reject;
153 
154  }
155 
157 
158 } // end namespace ifdb
Filter events based on their run/event numbers.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
std::unique_ptr< BeamFolder > fFolder
bool convertNovaTimeToUnixTime(uint64_t const &inputNovaTime, struct timespec &outputUnixTime)
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
Far Detector at Ash River, MN.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Reject triggers that are too close to beam spills.
RemoveBeamSpills(const fhicl::ParameterSet &pset)
bool filter(art::Event &evt)
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
Float_t e
Definition: plot.C:35
Encapsulate the geometry of one entire detector (near, far, ndos)
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11