test_fft.cxx
Go to the documentation of this file.
2 #include "WireCellUtil/Units.h"
4 
5 #include "TH1F.h"
6 #include "TLine.h"
7 #include "TText.h"
8 #include "TGraph.h"
9 
10 #include "MultiPdf.h"
11 
12 #include <iostream>
13 #include <algorithm>
14 #include <chrono>
15 
16 using namespace std;
17 using namespace WireCell;
18 using namespace WireCell::Test;
19 
20 // The preferred display units for gain.
21 const double GUnit = units::mV/units::fC;
22 
23 
25  Waveform::realseq_t& res,
26  const std::string& title,
27  const Binning& tbins)
28 {
31 
32  TH1F h_wave("response",title.c_str(),
33  tbins.nbins(), tbins.min()/units::us, tbins.max()/units::us);
34  TH1F h_wave2("response2",title.c_str(),
35  tbins.nbins(), tbins.min()/units::us, tbins.max()/units::us);
36  cout<<"tbins.min()"<<tbins.min()<<endl;
37  cout<<"tbins.max()"<<tbins.max()<<endl;
38  h_wave2.SetLineColor(2);
39 
40  h_wave.SetXTitle("Time (microsecond)");
41  h_wave.SetYTitle("Gain (mV/fC)");
42  h_wave.GetXaxis()->SetRangeUser(0,10.0);
43 
44  const int nticks = tbins.nbins();
45  for (int ind=0; ind<nticks; ++ind) {
46  h_wave.Fill(tbins.center(ind)/units::us, res[ind]/GUnit);
47  h_wave2.Fill(tbins.center(ind)/units::us, res2[ind]/GUnit);
48  }
49 
50  cerr << nticks << " " << spec.size() << endl;
51 
52  const double tick = tbins.binsize();
53  const Binning fbins(nticks, 0, 1/tick);
54 
55  TH1F h_mag("mag","Magnitude of Fourier transform of response",
56  fbins.nbins(), fbins.min()/units::megahertz, fbins.max()/units::megahertz);
57  h_mag.SetYTitle("Amplitude [mV/fC]");
58  h_mag.SetXTitle("MHz");
59 
60  TH1F h_phi("phi","Phase of Fourier transform of response",
61  fbins.nbins(), fbins.min()/units::megahertz, fbins.max()/units::megahertz);
62  h_phi.SetYTitle("Phase [radian]");
63  h_phi.SetXTitle("MHz");
64 
65  for (int ind=0; ind<nticks; ++ind) {
66  auto c = spec[ind];
67  const double freq = fbins.center(ind);
68  h_mag.Fill(freq/units::megahertz, std::abs(c)/GUnit);
69  h_phi.Fill(freq/units::megahertz, std::arg(c));
70  }
71  //h_mag.GetXaxis()->SetRangeUser(0,2.0);
72  //h_phi.GetXaxis()->SetRangeUser(0,2.0);
73 
74  pdf.canvas.Clear();
75  pdf.canvas.Divide(2,1);
76 
77  auto pad = pdf.canvas.cd(1);
78  pad->SetGridx();
79  pad->SetGridy();
80  h_wave.Draw("hist");
81  h_wave2.Draw("hist,same");
82 
83  auto spad = pdf.canvas.cd(2);
84  spad->Divide(1,2);
85  pad = spad->cd(1);
86  pad->SetGridx();
87  pad->SetGridy();
88  h_mag.Draw("hist");
89 
90  TLine ndb;
91  TText ndbtxt;
92  double maxmag = h_mag.GetMaximum();
93  double db3 = maxmag / 1.414;
94  double db6 = maxmag / 1.995;
95  double db10 = maxmag / 3.162;
96  double db20 = maxmag / 10;
97  ndb.DrawLine(0,db3, 1.0, db3);
98  ndb.DrawLine(0,db6, 1.0, db6);
99  ndb.DrawLine(0,db10, 1.0, db10);
100  ndb.DrawLine(0,db20, 1.0, db20);
101  ndbtxt.DrawText(0.5, db3, "-3dB");
102  ndbtxt.DrawText(0.5, db6, "-6dB");
103  ndbtxt.DrawText(0.5, db10, "-10dB");
104  ndbtxt.DrawText(0.5, db20, "-20dB");
105 
106  pad = spad->cd(2);
107  pad->SetGridx();
108  pad->SetGridy();
109  h_phi.Draw("hist");
110  pdf();
111 
112 }
113 
114 int main(int argc, char* argv[])
115 {
116  const std::vector<double> gains = {7.8*GUnit, 14.0*GUnit};
117  const std::vector<double> shapings = {1.0*units::us, 2.0*units::us};
118 
119  const double maxtime = 100.0*units::us;
120  const double tick = 0.5*units::us;
121  const Binning tbins(maxtime/tick, 0, maxtime);
122  const int nticks = tbins.nbins();
123  cerr << "Using nticks=" << nticks << endl;
124 
125  MultiPdf pdf(argv[0]);
126 
127  for (size_t ind=0; ind<gains.size(); ++ind) {
128  Response::ColdElec ce(gains[ind], shapings[ind]);
129  Waveform::realseq_t res = ce.generate(tbins);
130 
131  const double tshape_us = shapings[ind]/units::us;
132  auto tit = Form("Cold Electronics Response at %.0fus peaking",
133  tshape_us);
134  draw_time_freq(pdf, res, tit, tbins);
135  }
136 
137 
138  // Look at RC filter
139 
140  const double tconst = 1000.0*units::ms;
141 
142  {
143  Response::SimpleRC rc(tconst);
144  Waveform::realseq_t res = rc.generate(tbins);
145 
146  auto tit = "RC Response at 1ms time constant";
147  draw_time_freq(pdf, res, tit, tbins);
148  }
149  {
150  Binning shifted(tbins.nbins(), tbins.min()+tick, tbins.max()+tick);
151 
152  Response::SimpleRC rc(tconst);
153  Waveform::realseq_t res = rc.generate(shifted);
154 
155  auto tit = "RC Response at 1ms time constant (suppress delta)";
156  draw_time_freq(pdf, res, tit, tbins);
157  }
158 
159  // Look at SysResp (Gaussian smear)
160  {
162  Waveform::realseq_t res = gaus.generate(tbins);
163  auto tit = "Response Gaussian smear by default";
164  draw_time_freq(pdf, res, tit, tbins);
165  }
166  {
167  double mag = 1.0;
168  double smear = 2.0*units::us;
169  double nsigma = 5;
170  Binning ttt(nsigma*2*smear/tick, -nsigma*smear, nsigma*smear);
171  //Binning ttt(nsigma*2*smear/tick, 0, 2*nsigma*smear);
172  Response::SysResp gaus(tick, mag, smear);
173  Waveform::realseq_t res = gaus.generate(ttt);
174  auto tit = "Response Gaussian 2 us smear";
175  draw_time_freq(pdf, res, tit, ttt);
176  }
177 
178 
179  // do timing tests
180  {
181 
182 
183  TGraph* timings[4] = {
184  new TGraph, new TGraph, new TGraph, new TGraph
185  };
186 
187  // Some popular choices with powers-of-two sprinkled in
188  std::vector<int> nsampleslist{128, 256,
189  400,480, // protoDUNE U/V and W channels per plane
190  512,
191  800, // protoDUNE, sum of U or V channels for both faces
192  960, // protoDUNE, sum of W channels (or wires) for both faces
193  1024,
194  1148, // N wires in U/V plane for protodune
195  2048,
196  2400, // number of channels in U or V in microboone
197  2560, // DUNE, total APA channels
198  3456, // number of channels in microboone's W
199  4096,
200  6000, // one choice of nticks for protoDUNE
201  8192,
202  8256, // total microboone channels
203  9592,9594,9595,9600, // various microboone readout lengths
204  10000, // 5 ms at 2MHz readout
205  10240,
206  16384};
207  const int ntries = 1000;
208  for (auto nsamps : nsampleslist) {
209  Response::ColdElec ce(gains[1], shapings[1]);
210  const Binning bins(nsamps, 0, maxtime);
211  Waveform::realseq_t res = ce.generate(bins);
212  Waveform::compseq_t spec;
213 
214  double fwd_time = 0.0;
215  for (int itry=0; itry<ntries; ++itry) {
217  spec = Waveform::dft(res);
219  fwd_time += std::chrono::duration_cast<std::chrono::nanoseconds>(t2-t1).count();
220  }
221  fwd_time /= ntries;
222 
223  double rev_time = 0.0;
224  for (int itry=0; itry<ntries; ++itry) {
226  res = Waveform::idft(spec);
228  rev_time += std::chrono::duration_cast<std::chrono::nanoseconds>(t2-t1).count();
229  }
230  rev_time /= ntries;
231 
232  cerr << "DFT nsamples=" << nsamps
233  << "\n\tforward: " << fwd_time/1000.0 << " us, " << fwd_time/nsamps/1000.0 << " us/sample"
234  << "\n\treverse: " << rev_time/1000.0 << " us, " << rev_time/nsamps/1000.0 << " us/sample"
235  << "\n\taverage: " << 0.5*(fwd_time+rev_time)/1000.0 << " us"
236  << endl;
237 
238  cout << "timing " << nsamps << " "
239  << fwd_time/1000.0 << " "
240  << rev_time/1000.0 << "\n";
241 
242  timings[0]->SetPoint(timings[0]->GetN(), nsamps, fwd_time);
243  timings[1]->SetPoint(timings[1]->GetN(), nsamps, fwd_time/nsamps);
244  timings[2]->SetPoint(timings[2]->GetN(), nsamps, rev_time);
245  timings[3]->SetPoint(timings[3]->GetN(), nsamps, rev_time/nsamps);
246  }
247 
248  pdf.canvas.Clear();
249  pdf.canvas.Divide(1,2);
250 
251  auto text = new TText;
252  {
253  auto pad = pdf.canvas.cd(1);
254  pad->SetGridx();
255  pad->SetGridy();
256  pad->SetLogx();
257  auto graph = timings[0];
258  auto frame = graph->GetHistogram();
259  frame->SetTitle("Fwd/rev DFT timing (absolute)");
260  frame->GetXaxis()->SetTitle("number of samples");
261  frame->GetYaxis()->SetTitle("time (ns)");
262  timings[0]->Draw("AL");
263  timings[2]->Draw("L");
264  for (int ind=0; ind<graph->GetN(); ++ind) {
265  auto x = graph->GetX()[ind];
266  auto y = graph->GetY()[ind];
267  text->DrawText(x,y,Form("%.0f", x));
268  }
269  }
270 
271  {
272  auto pad = pdf.canvas.cd(2);
273  pad->SetGridx();
274  pad->SetGridy();
275  pad->SetLogx();
276  auto frame = timings[1]->GetHistogram();
277  frame->SetTitle("Fwd/rev DFT timing (relative)");
278  frame->GetXaxis()->SetTitle("number of samples");
279  frame->GetYaxis()->SetTitle("time per sample (ns/samp)");
280  timings[1]->Draw("AL");
281  timings[3]->Draw("L");
282  }
283  pdf();
284 
285  }
286 
287  return 0;
288 }
Sequence< real_t > realseq_t
Definition: Waveform.h:31
double center(int ind) const
Definition: Binning.h:86
std::string string
Definition: nybbler.cc:12
internal::named_arg< T, char > arg(string_view name, const T &arg)
Definition: core.h:1391
def graph(desc, maker=maker)
Definition: apa.py:294
STL namespace.
A functional object caching gain and shape.
Definition: Response.h:165
const double tick
const int nticks
double binsize() const
Definition: Binning.h:73
double max() const
Definition: Binning.h:52
Binning tbins(nticks, t0, t0+readout_time)
double y
void draw_time_freq(MultiPdf &pdf, Waveform::realseq_t &res, const std::string &title, const Binning &tbins)
Definition: test_fft.cxx:24
T abs(T value)
T gaus(T x, T amplitude, T mean, T sigma)
Tests GausssianFit object with a known input.
float fC
Definition: units.py:113
WireCell::Waveform::realseq_t generate(const WireCell::Waveform::Domain &domain, int nsamples)
FIXME: eradicate Domain in favor of Binning.
Definition: Response.cxx:303
double min() const
Definition: Binning.h:47
megahertz_as<> megahertz
Type of frequency stored in megahertz, in double precision.
Definition: frequency.h:101
const double GUnit
Definition: test_fft.cxx:21
Definition: Main.h:22
static const double us
Definition: Units.h:101
static const double ms
Definition: Units.h:100
int main(int argc, char *argv[])
Definition: test_fft.cxx:114
list x
Definition: train.py:276
array_xxf idft(const array_xxc &arr)
Definition: Array.cxx:96
array_xxc dft(const array_xxf &arr)
Definition: Array.cxx:15
int nbins() const
Definition: Binning.h:42
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34
QTextStream & endl(QTextStream &s)
nanosecond nanoseconds
Alias for common language habits.
Definition: spacetime.h:134