Functions
test_convo.cxx File Reference
#include "WireCellUtil/Response.h"
#include "WireCellUtil/Waveform.h"
#include "WireCellUtil/Units.h"
#include "WireCellUtil/ExecMon.h"
#include "WireCellUtil/Testing.h"
#include "WireCellUtil/Binning.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TLine.h"
#include "TStyle.h"
#include "TH2F.h"
#include <vector>
#include <iostream>
#include <map>

Go to the source code of this file.

Functions

std::vector< TH1F * > plot_wave (TCanvas &canvas, int padnum, std::string name, std::pair< double, double > bounds, const Binning &bins, const Waveform::realseq_t &tdomain, const Waveform::compseq_t &fdomain)
 
int main (int argc, char *argv[])
 

Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 72 of file test_convo.cxx.

73 {
74  if (argc < 2) {
75  std::cerr << "This test requires an Wire Cell Field Response input file." << std::endl;
76  return 0;
77  }
78 
79  const char* uvw="UVW";
80 
81  // time binning
82  const int nticks = 10000;
83  const double tick = 0.5*units::us;
84  Binning bins(nticks, 0, nticks*tick);
85 
86  // Get one field response (the one nearest the WOI)
87  auto fr = Response::Schema::load(argv[1]);
88  cerr << "Drift speed is: " << fr.speed/(units::mm/units::us) << " mm/us\n";
89 
90  TCanvas canvas("h","h",900,1200);
91  canvas.Print(Form("%s.pdf[", argv[0]), "pdf");
92 
93  const int nplanes = 3;
94  for (int iplane=0; iplane<nplanes; ++iplane) {
95 
96  auto plane_response = fr.planes[iplane];
97 
98  const int npaths = plane_response.paths.size();
99 
100  for (int ipath=0; ipath<npaths; ++ipath) {
101  auto& path_response = plane_response.paths[ipath];
102 
103  WireCell::Waveform::realseq_t raw_response = path_response.current;
104 
105  const int rawresp_size = raw_response.size();
106  Assert(rawresp_size);
107  const double rawresp_min = fr.tstart;
108  const double rawresp_tick = fr.period;
109  const double rawresp_max = rawresp_min + rawresp_size*rawresp_tick;
110  Binning rawresp_bins(rawresp_size, rawresp_min, rawresp_max);
111 
112  // match response sampling to digi and zero-pad
114  for (int ind=0; ind<rawresp_size; ++ind) {
115  const double time = rawresp_bins.center(ind);
116  const int bin = bins.bin(time);
117  response[bin] += raw_response[ind];
118  }
119 
120  // Make up a Gaussian charge distribution
121  const double charge_const = 1000.0;
122  const double charge_time = 1*units::ms;
123  const double charge_sigma = 3*units::us;
124 
125  WireCell::Waveform::realseq_t electrons(nticks, 0.0);
126  for (int ind=0; ind<nticks; ++ind) {
127  const double t = bins.center(ind);
128  const double rel = (t-charge_time)/charge_sigma;
129  const double val = charge_const * exp(-0.5*rel*rel);
130  electrons[ind] = val;
131  }
132 
133  // frequency space
134  Waveform::compseq_t charge_spectrum = Waveform::dft(electrons);
135  Waveform::compseq_t raw_response_spectrum = Waveform::dft(raw_response);
136  Waveform::compseq_t response_spectrum = Waveform::dft(response);
137 
138  // convolve
139  Waveform::compseq_t conv_spectrum(nticks, Waveform::complex_t(0.0,0.0));
140  for (int ind=0; ind < nticks; ++ind) {
141  conv_spectrum[ind] = response_spectrum[ind]*charge_spectrum[ind];
142  }
143  Waveform::realseq_t conv = Waveform::idft(conv_spectrum);
144  for (int ind=0; ind < nticks; ++ind) {
145  conv[ind] /= nticks;
146  }
147 
148  canvas.Clear();
149  canvas.Divide(3,4);
150 
151  std::string extra = Form(" %c-%.1f", uvw[iplane], path_response.pitchpos);
152 
153  plot_wave(canvas, 1, "Response"+extra, std::make_pair(0,.1),
154  rawresp_bins, raw_response, raw_response_spectrum);
155  plot_wave(canvas, 4, "Response"+extra, std::make_pair(0,.1),
156  bins, response, response_spectrum);
157  plot_wave(canvas, 7, "Electrons"+extra, std::make_pair(.9,1.1),
158  bins, electrons, charge_spectrum);
159  plot_wave(canvas, 10, "Signal"+extra, std::make_pair(.9,1.1),
160  bins, conv, conv_spectrum);
161 
162  canvas.Print(Form("%s.pdf", argv[0]), "pdf");
163 
164  } // paths
165  } // planes
166 
167  canvas.Print(Form("%s.pdf]", argv[0]), "pdf");
168  return 0;
169 }
Sequence< real_t > realseq_t
Definition: Waveform.h:31
std::string string
Definition: nybbler.cc:12
const double tick
const int nticks
#define Assert
Definition: Testing.h:7
std::vector< TH1F * > plot_wave(TCanvas &canvas, int padnum, std::string name, std::pair< double, double > bounds, const Binning &bins, const Waveform::realseq_t &tdomain, const Waveform::compseq_t &fdomain)
Definition: test_convo.cxx:23
static const double mm
Definition: Units.h:73
static const double us
Definition: Units.h:101
static const double ms
Definition: Units.h:100
QTextStream & bin(QTextStream &s)
array_xxf idft(const array_xxc &arr)
Definition: Array.cxx:96
array_xxc dft(const array_xxf &arr)
Definition: Array.cxx:15
std::complex< float > complex_t
The type for the spectrum in each bin.
Definition: Waveform.h:21
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34
QTextStream & endl(QTextStream &s)
def load(filename, jpath="depos")
Definition: depos.py:34
std::vector<TH1F*> plot_wave ( TCanvas &  canvas,
int  padnum,
std::string  name,
std::pair< double, double >  bounds,
const Binning bins,
const Waveform::realseq_t tdomain,
const Waveform::compseq_t fdomain 
)

Definition at line 23 of file test_convo.cxx.

28 {
29  const int nbins = bins.nbins();
30  const double vmax = (1.0/(bins.binsize()/units::second))*1e-6;
31  const double vbin = vmax/nbins;
32 
33  TH1F* ht = new TH1F("ht", Form("%s - waveform", name.c_str()),
34  nbins, bins.min()/units::ms, bins.max()/units::ms);
35  ht->SetXTitle(Form("time (ms, %d bins)", nbins));
36  ht->SetDirectory(0);
37  ht->GetXaxis()->SetRangeUser(bounds.first, bounds.second);
38 
39  TH1F* hm = new TH1F("hm", Form("%s - magnitude", name.c_str()),
40  nbins, 0, vmax);
41  hm->SetXTitle(Form("frequency (MHz, %d bins)", nbins));
42  hm->SetDirectory(0);
43 
44  TH1F* hp = new TH1F("hp", Form("%s - phase", name.c_str()),
45  nbins, 0, vmax);
46  hp->SetXTitle(Form("frequency (MHz, %d bins)", nbins));
47  hp->SetDirectory(0);
48 
49  for (int ind=0; ind<bins.nbins(); ++ind) {
50  const double tms = bins.center(ind);
51  const double fMHz = (ind+0.5)*vbin;
52  ht->Fill(tms/units::ms, tdomain[ind]);
53  hm->Fill(fMHz, std::abs(fdomain[ind]));
54  hp->Fill(fMHz, std::arg(fdomain[ind]));
55  }
56 
57  TVirtualPad* gpad=0;
58 
59  gpad = canvas.cd(padnum++);
60  ht->Draw("HIST");
61 
62  gpad = canvas.cd(padnum++);
63  gpad->SetLogy(1);
64  hm->Draw("HIST");
65 
66  gpad = canvas.cd(padnum++);
67  hp->Draw("HIST");
68 
69  return std::vector<TH1F*>{ht,hm,hp};
70 }
static QCString name
Definition: declinfo.cpp:673
double center(int ind) const
Definition: Binning.h:86
internal::named_arg< T, char > arg(string_view name, const T &arg)
Definition: core.h:1391
double binsize() const
Definition: Binning.h:73
double max() const
Definition: Binning.h:52
T abs(T value)
const double e
double min() const
Definition: Binning.h:47
BoundingBox bounds(int x, int y, int w, int h)
Definition: main.cpp:37
static const double ms
Definition: Units.h:100
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:80
int nbins() const
Definition: Binning.h:42