|
| 1 | +#include <iostream> |
| 2 | +#include <stdlib.h> |
| 3 | +#include <string> |
| 4 | +#include <vector> |
| 5 | +#include <numeric> |
| 6 | +#include <getopt.h> |
| 7 | +#include <chrono> |
| 8 | +#include <float.h> |
| 9 | + |
| 10 | + |
| 11 | +//some ROOT includes |
| 12 | +#include "TInterpreter.h" |
| 13 | +#include "TROOT.h" |
| 14 | +#include "TH1F.h" |
| 15 | +#include "TH2D.h" |
| 16 | +#include "TCanvas.h" |
| 17 | +#include "TTree.h" |
| 18 | +#include "TFile.h" |
| 19 | +#include "TStyle.h" |
| 20 | +#include "TSystem.h" |
| 21 | +#include "TGraph.h" |
| 22 | +#include "TFFTReal.h" |
| 23 | +#include "TTimeStamp.h" |
| 24 | +#include "TLatex.h" |
| 25 | +#include "TLine.h" |
| 26 | + |
| 27 | + |
| 28 | +// art includes |
| 29 | +#include "canvas/Utilities/InputTag.h" |
| 30 | +#include "canvas/Persistency/Common/FindOneP.h" |
| 31 | +#include "canvas/Persistency/Common/FindManyP.h" |
| 32 | +#include "canvas/Persistency/Common/Ptr.h" |
| 33 | +#include "art/Framework/Core/EDFilter.h" |
| 34 | +#include "art/Framework/Core/ModuleMacros.h" |
| 35 | +#include "art/Framework/Principal/Event.h" |
| 36 | +#include "art/Framework/Principal/Handle.h" |
| 37 | +#include "art/Framework/Principal/Run.h" |
| 38 | +#include "art/Framework/Principal/SubRun.h" |
| 39 | +#include "fhiclcpp/ParameterSet.h" |
| 40 | +#include "messagefacility/MessageLogger/MessageLogger.h" |
| 41 | +#include "lardataobj/RecoBase/Hit.h" |
| 42 | +#include "fhiclcpp/ParameterSet.h" |
| 43 | +#include "lardataobj/RawData/RawDigit.h" |
| 44 | +#include "lardataobj/RawData/RDTimeStamp.h" |
| 45 | +#include "lardataobj/RawData/raw.h" |
| 46 | +#include "larcore/Geometry/Geometry.h" |
| 47 | +#include "canvas/Persistency/Provenance/Timestamp.h" |
| 48 | + |
| 49 | +#include "art_root_io/TFileService.h" |
| 50 | + |
| 51 | +#include "sbndcode/ChannelMaps/TPC/TPCChannelMapService.h" |
| 52 | + |
| 53 | + |
| 54 | +namespace filt{ |
| 55 | + |
| 56 | + class FilterChoppy : public art::EDFilter { |
| 57 | + public: |
| 58 | + |
| 59 | + explicit FilterChoppy(fhicl::ParameterSet const & pset); |
| 60 | + virtual bool filter(art::Event& evt) override; |
| 61 | + void reconfigure(fhicl::ParameterSet const& pset); |
| 62 | + virtual void beginJob() override; |
| 63 | + |
| 64 | + std::vector<art::Ptr<raw::RDTimeStamp>> raw_timestamps_handle; |
| 65 | + |
| 66 | + private: |
| 67 | + |
| 68 | + int run; |
| 69 | + int subrun; |
| 70 | + int event; |
| 71 | + long delta_tmaxmin; |
| 72 | + std::vector<long> ts_vec; |
| 73 | + std::vector<long> ts_hist_vec; |
| 74 | + int npeak; |
| 75 | + |
| 76 | + }; |
| 77 | + |
| 78 | + |
| 79 | + FilterChoppy::FilterChoppy(fhicl::ParameterSet const & pset) |
| 80 | + : EDFilter(pset) |
| 81 | + { |
| 82 | + this->reconfigure(pset); |
| 83 | + } |
| 84 | + |
| 85 | + void FilterChoppy::reconfigure(fhicl::ParameterSet const& pset){ |
| 86 | + } |
| 87 | + |
| 88 | + |
| 89 | + bool FilterChoppy::filter(art::Event & evt){ |
| 90 | + |
| 91 | + raw_timestamps_handle.clear(); |
| 92 | + run = evt.id().run(); |
| 93 | + subrun = evt.id().subRun(); |
| 94 | + event = evt.id().event(); |
| 95 | + |
| 96 | + ts_vec.clear(); |
| 97 | + ts_hist_vec.clear(); |
| 98 | + npeak = 0; |
| 99 | + delta_tmaxmin = 0; |
| 100 | + |
| 101 | + // get the timestamps |
| 102 | + art::Handle<std::vector<raw::RDTimeStamp>> timestamp_handle; |
| 103 | + evt.getByLabel("daq", timestamp_handle); |
| 104 | + |
| 105 | + // exit if the data isn't present |
| 106 | + if (!timestamp_handle.isValid()) { |
| 107 | + std::cerr << "Error: missing timestamps with producer (" << "daq" << ")" << std::endl; |
| 108 | + return false; |
| 109 | + } |
| 110 | + art::fill_ptr_vector(raw_timestamps_handle, timestamp_handle); |
| 111 | + |
| 112 | + // collect TS |
| 113 | + unsigned timeindex = 0; |
| 114 | + for (auto const& timestamps: raw_timestamps_handle) { |
| 115 | + long this_ts = timestamps->GetTimeStamp(); |
| 116 | + ts_vec.push_back(this_ts); |
| 117 | + timeindex++; |
| 118 | + } |
| 119 | + // find min, max TS values |
| 120 | + long ts_min = ts_vec[0]; |
| 121 | + long ts_max = ts_vec[0]; |
| 122 | + for (long ts : ts_vec) { |
| 123 | + if (ts < ts_min) { |
| 124 | + ts_min = ts; |
| 125 | + } |
| 126 | + if (ts > ts_max) { |
| 127 | + ts_max = ts; |
| 128 | + } |
| 129 | + } |
| 130 | + delta_tmaxmin = std::abs(ts_max - ts_min); |
| 131 | + |
| 132 | + // if (delta_tmaxmin > 100000){ |
| 133 | + // std::cout << "Event " << event << " has a delta_tmaxmin of " << delta_tmaxmin << std::endl; |
| 134 | + // } |
| 135 | + return delta_tmaxmin < 100000; |
| 136 | + |
| 137 | + |
| 138 | + // // Identify choppy events by looking at the # of peaks in the TS distribution |
| 139 | + // make binned hist |
| 140 | + // float nbin = 100.; |
| 141 | + // float bin_width = (ts_max - ts_min)/nbin; |
| 142 | + // for (int i = 0; i < 100; ++i) { |
| 143 | + // float bin_low = ts_min + i*bin_width; |
| 144 | + // float bin_high = ts_min + (i+1)*bin_width; |
| 145 | + // int bin_count = 0; |
| 146 | + // for (long ts : ts_vec) { |
| 147 | + // if ((bin_low < ts) & (ts <= bin_high)) { |
| 148 | + // bin_count += 1; |
| 149 | + // } |
| 150 | + // } |
| 151 | + // ts_hist_vec.push_back(bin_count); |
| 152 | + // } |
| 153 | + // // find peak bins |
| 154 | + // for (int bc : ts_hist_vec) { |
| 155 | + // //TODO: make threshold fcl parameter |
| 156 | + // if (bc > 1000) { |
| 157 | + // npeak += 1; |
| 158 | + // } |
| 159 | + // } |
| 160 | + |
| 161 | + // // If the energy deposit within the beam time is greater than some limit then trigger the event |
| 162 | + // if (npeak > 1){ |
| 163 | + // std::cout << "Event " << event << " has " << npeak << " peaks" << std::endl; |
| 164 | + // } |
| 165 | + // return npeak < 2; |
| 166 | + |
| 167 | + } |
| 168 | + |
| 169 | + void FilterChoppy::beginJob() { |
| 170 | + } |
| 171 | + |
| 172 | + DEFINE_ART_MODULE(FilterChoppy) |
| 173 | + |
| 174 | +} |
0 commit comments