EventWaveformDump_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \file EventWaveformDump_module.cc
3 /// \module analyzer
4 /// \brief Dump waveforms for a range of events
5 /// \author Aidan Medcalf (University of Dallas) <amedcalf@fnal.gov>
6 /// \date May 2019
7 ////////////////////////////////////////////////////////////////////////////
8 
9 // framework
13 #include "fhiclcpp/ParameterSet.h"
21 
22 // nova
23 #include "RawData/RawBeamline.h"
26 #include "BeamlineRecoBase/ToF.h"
29 
30 // ROOT
31 #include "TH1D.h"
32 #include "TH2D.h"
33 #include "TCanvas.h"
34 #include "TLegend.h"
35 #include "TFile.h"
36 #include "TGraph.h"
37 #include "TMultiGraph.h"
38 #include "THStack.h"
39 #include "TLine.h"
40 #include "TText.h"
41 
42 // stl
43 #include <iostream>
44 
45 // root
46 
47 // -----------------------------------------------------------------------
48 namespace beamlinereco {
49 
51 
52  public:
53 
55 
56  void reconfigure(const fhicl::ParameterSet& pset);
57  void beginJob();
58  void analyze(const art::Event& e);
59 
60  private:
66  std::vector<std::pair<unsigned int, std::string>> fDumpChannels;
69  bool fDrawToFs;
70 
71  std::map<unsigned int, std::string> dumpChannelMap;
72 
73  // "Persistence plot" of all waveforms overlayed, per channel
74  std::map<unsigned int, TH2D*> hPersistence;
75 
76  // services
80 
81  unsigned int onlineChannel(ChannelID id);
82  };
83 
85 
86 }
87 
89  if (id.System == 2) {
90  switch (id.Detector) {
91  case 0: return 8 + id.Channel;
92  case 1: return 12 + id.Channel;
93  case 2: return 4 + id.Channel;
94  default: return 99999;
95  }
96  }
97  return 99990;
98 }
99 
100 // -----------------------------------------------------------------------
102  this->reconfigure(pset);
103 }
104 
105 // -----------------------------------------------------------------------
107  fRawToFDigitLabel = pset.get<art::InputTag>("RawToFDigitLabel");
108  fRawCherenkovDigitLabel = pset.get<art::InputTag>("RawCherenkovDigitLabel");
109  fRecoToFDigitLabel = pset.get<art::InputTag>("RecoToFDigitLabel");
110  fRecoCherenkovDigitLabel = pset.get<art::InputTag>("RecoCherenkovDigitLabel");
111  fRecoToFLabel = pset.get<art::InputTag>("RecoToFLabel");
112  fDumpChannels = pset.get<std::vector<std::pair<unsigned int, std::string>>>("DumpChannels");
113  fSaveWaveforms = pset.get<bool>("SaveWaveforms");
114  fDrawMarkers = pset.get<bool>("DrawMarkers");
115  fDrawToFs = pset.get<bool>("DrawToFs");
116 
117  dumpChannelMap.clear();
118  for (auto ch : fDumpChannels) dumpChannelMap[ch.first] = ch.second;
119 }
120 
121 // -----------------------------------------------------------------------
123  art::TFileDirectory dir = tfs->mkdir("Persistence");
124  for (auto ch : fDumpChannels) {
125  hPersistence[ch.first] = dir.make<TH2D>(Form("hPersistence_channel%d",ch.first),
126  Form("Persistence (Channel %d);Time [ticks];ADC;",ch.first),
127  1024, 0, 1024,
128  4096, 0, 4096);
129  }
130 }
131 
132 // -----------------------------------------------------------------------
134  //std::cout << "+++++++++++++++Run " << e.run() << " Subrun " << e.subRun() << " Event " << e.event() << std::endl;
135 
136  std::vector<art::Ptr<rawdata::RawBeamlineDigit>> digits;
137  std::vector<art::Ptr<rawdata::RawBeamlineDigit>> tofDigits;
138  std::vector<art::Ptr<rawdata::RawBeamlineDigit>> cherenkovDigits;
139 
140  std::vector<art::Ptr<brb::BeamlineDigit>> recoDigits;
141  std::vector<art::Ptr<brb::BeamlineDigit>> recoToFDigits;
142  std::vector<art::Ptr<brb::BeamlineDigit>> recoCherenkovDigits;
143 
144  // Get raw tof digits
146  if (e.getByLabel(fRawToFDigitLabel, tofDigitHandle))
147  art::fill_ptr_vector(tofDigits, tofDigitHandle);
148  // Get raw cherenkov digits
150  if (e.getByLabel(fRawCherenkovDigitLabel, cherenkovDigitHandle))
151  art::fill_ptr_vector(cherenkovDigits, cherenkovDigitHandle);
152  // Combine raw vectors
153  digits.reserve(tofDigits.size() + cherenkovDigits.size());
154  digits.insert(digits.end(), tofDigits.begin(), tofDigits.end());
155  digits.insert(digits.end(), cherenkovDigits.begin(), cherenkovDigits.end());
156 
157  // Get reco tof digits
159  if (e.getByLabel(fRecoToFDigitLabel, recoToFDigitHandle))
160  art::fill_ptr_vector(recoToFDigits, recoToFDigitHandle);
161  // Get reco cherenkov digits
162  art::Handle<std::vector<brb::BeamlineDigit>> recoCherenkovDigitHandle;
163  if (e.getByLabel(fRecoCherenkovDigitLabel, recoCherenkovDigitHandle))
164  art::fill_ptr_vector(recoCherenkovDigits, recoCherenkovDigitHandle);
165  // Combine reco vectors
166  recoDigits.reserve(recoToFDigits.size() + recoCherenkovDigits.size());
167  recoDigits.insert(recoDigits.end(), recoToFDigits.begin(), recoToFDigits.end());
168  recoDigits.insert(recoDigits.end(), recoCherenkovDigits.begin(), recoCherenkovDigits.end());
169 
170  // Get reco tofs
172  std::vector<art::Ptr<brb::ToF>> tofs;
173  if (e.getByLabel(fRecoToFLabel, tofHandle))
174  art::fill_ptr_vector(tofs, tofHandle);
175 
176  // If we have no digits, exit early
177  if (digits.size() < 1) return;
178 
179  bool haveChannels = false;
180  for (auto d : digits) {
181  // if channel is in the dump list
182  if (dumpChannelMap.count(onlineChannel(d->ChannelID()))) {
183  haveChannels = true;
184  break;
185  }
186  }
187  if (!haveChannels) return;
188 
189  for (auto digit : digits) {
190  unsigned int signalchannel = onlineChannel(digit->ChannelID());
191  //std::cout << "Run " << e.run() << " Subrun " << e.subRun() << " Event " << e.event() << " OnlineDigitChannel " << signalchannel << " " << digit->ChannelID() << std::endl;
192  if (dumpChannelMap.count(signalchannel) == 0) continue;
193  //std::cout << "Channel: " << signalchannel << " about to fill persistence" << std::endl;
194  for (unsigned int i=0; i<digit->NADC(); i++) {
195  hPersistence[signalchannel]->Fill(i, digit->ADC(i));
196  }
197  //std::cout << "Filled" << std::endl;
198  }
199 
200  art::TFileDirectory dir = tfs->mkdir(Form("Run%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()));
201  art::TFileDirectory hdir = dir.mkdir("Histograms");
202 
203  TLegend *leg = new TLegend(0.70,0.20,0.88,0.60);
204 
205  std::vector<TH1D*> waveforms;
206  std::vector<TGraph*> markers;
207 
208  //std::cout << "------------------------------"
209  //<< "Run " << e.run() << " Subrun " << e.subRun() << " Event " << e.event() << std::endl;
210  //for (auto it : dumpChannelMap) { std::cout << "\t\t[" << it.first << ":" << it.second << "]" << std::endl; }
211  // For each digit
212  for (auto digit : digits) {
213  unsigned int signalchannel = onlineChannel(digit->ChannelID());
214  //std::cout << signalchannel << " " << dumpChannelMap.count(signalchannel) << std::endl;
215  if (dumpChannelMap.count(signalchannel) == 0) continue;
216 
217  // Create histogram for waveform
218  TH1D *waveform = hdir.make<TH1D>(Form("hr%ds%de%d_ch%d",e.run(),e.subRun(),e.event(),signalchannel), ";Time [ticks];ADC;", digit->NADC(), 0, digit->NADC());
219  //std::cout << "channel: " << signalchannel << "\t\t" << digit->ChannelID() << std::endl;
220  //std::cout << "\t\t\tDigit has " << digit->NADC() << std::endl;
221  //std::cout << signalchannel << " [";
222  // Fill
223  for (unsigned int i=0; i<digit->NADC(); i++) {
224  //std::cout << digit->ADC(i) << ", ";
225  waveform->Fill(i, digit->ADC(i));
226  }
227  //std::cout << "]" << std::endl;
228 
229  //std::cout << "Filled waveform" << std::endl;
230 
231  // Draw
232  waveform->GetYaxis()->SetRangeUser(0,5000);
233  waveform->SetLineWidth(2);
234  leg->AddEntry(waveform, Form("Channel %d %s", signalchannel, dumpChannelMap[signalchannel].c_str()), "l");
235  waveforms.push_back(waveform);
236 
237  if (fDrawMarkers) {
238  std::vector<art::Ptr<brb::BeamlineDigit>> rds;
239  for (auto rd : recoDigits) { if (onlineChannel(rd->ChannelID()) == signalchannel) rds.push_back(rd); }
240  if (rds.size() > 0) {
241  TGraph *g = new TGraph(rds.size() * 2);
242  for (size_t i = 0; i<rds.size(); i++) {
243  double starttime = rds[i]->StartTimeInNanoSec() - fCalibration->TimeCorrection(rds[i]->ChannelID());
244  double peaktime = rds[i]->PeakTimeInNanoSec() - fCalibration->TimeCorrection(rds[i]->ChannelID());
245  //double starttime = rds[i]->StartTimeInNanoSec();
246  //double peaktime = rds[i]->PeakTimeInNanoSec();
247  unsigned int startidx = std::min((unsigned int)1023, (unsigned int)std::max(0.0, (starttime/0.4)));
248  unsigned int peakidx = std::min((unsigned int)1023, (unsigned int)std::max(0.0, (peaktime/0.4)));
249  //unsigned int startidx = (unsigned int)std::max(0.0, (starttime/0.4));
250  //unsigned int peakidx = (unsigned int)std::max(0.0, (peaktime/0.4));
251  //std::cout << "st=" << starttime << " si=" << startidx << " pt=" << peaktime << " pi=" << peakidx << std::endl;
252  double startADC = digit->ADC(startidx);
253  double peakADC = digit->ADC(peakidx);
254  g->SetPoint(2*i+0, starttime/0.4, startADC);
255  g->SetPoint(2*i+1, peaktime/0.4, peakADC);
256  }
257  g->SetMarkerColor(1);
258  g->SetMarkerStyle(3);
259  g->SetMarkerSize(2);
260  markers.push_back(g);
261  }
262  }
263  }
264 
265  // We manually write canvases to the TFileSystem file, which is bad form
266  TCanvas *canvas = dir.make<TCanvas>(Form("cRun%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()), "", 800, 600);
267  canvas->cd();
268 
269  THStack *stack = new THStack(Form("hsRun%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()), Form("hsRun%dSubRun%dEvent%d",e.run(),e.subRun(),e.event()));
270  for (auto w : waveforms) stack->Add(w, "hist");
271  stack->Draw("nostack hist PLC");
272 
273  if (fDrawMarkers) {
274  for (auto g : markers) g->Draw("same p");
275  }
276 
277  if (fDrawToFs) {
278  for (size_t i=0; i<tofs.size(); i++) {
279  auto tof = tofs[i];
280  TLine *l = new TLine(tof->Timestamps().first/0.4, 4900 - 200*i, tof->Timestamps().second/0.4, 4900 - 200*i);
281  l->SetLineColor(1);
282  l->SetLineWidth(2);
283  l->SetLineStyle(2);
284  l->Draw("same");
285  double tx = (tof->Timestamps().first + tof->Timestamps().second)/2.0/0.4;
286  double ty = 4900 - 200*i + 50;
287  TText *t = new TText(tx, ty, Form("tof = %f ns", tof->Time()));
288  t->SetTextAlign(21);
289  t->SetTextSize(0.02);
290  t->Draw();
291  }
292  }
293 
294  leg->Draw("same");
295  canvas->Write();
296 
297  delete leg;
298 }
T max(const caf::Proxy< T > &a, T b)
SubRunNumber_t subRun() const
Definition: Event.h:72
EventWaveformDump(const fhicl::ParameterSet &pset)
art::ServiceHandle< art::TFileService > tfs
art::ServiceHandle< beamlineutil::BeamlineChannelMap > fChannelMap
unsigned int onlineChannel(ChannelID id)
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
DEFINE_ART_MODULE(TestTMapFile)
list markers
Definition: styles.py:7
double TimeCorrection(ChannelID channel) const
Relative timing calibrations between individual digitizer channels.
Encapsulation of reconstructed digitizer &#39;hits&#39;. Used for ToF PMTs and SiPMs, and Cherenkov and Muon ...
std::map< unsigned int, TH2D * > hPersistence
T get(std::string const &key) const
Definition: ParameterSet.h:231
Encapsulation of reconstructed Time-of-Flight (ToF) information. Part of beamline reconstruction for ...
Float_t d
Definition: plot.C:236
EventNumber_t event() const
Definition: Event.h:67
Calibration service to provide specific externally-measured values to the relevant reconstruction/ana...
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
std::map< unsigned int, std::string > dumpChannelMap
art::ServiceHandle< beamlineutil::BeamlineCalibration > fCalibration
T * make(ARGS...args) const
void reconfigure(const fhicl::ParameterSet &pset)
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::vector< std::pair< unsigned int, std::string > > fDumpChannels
Raw data definitions for beamline data used in NOvA test beam experiment.
T min(const caf::Proxy< T > &a, T b)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Encapsulation of reconstructed PID information from detectors in the beamline (ToF, WCs, Cherenkov). Part of beamline reconstruction for NOvA test beam.
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
Float_t w
Definition: plot.C:20
Channel mapping service which may be used to interpret the channels which are read out by the various...