test_AdcDeconvoluteFFT.cxx
Go to the documentation of this file.
1 // test_AdcDeconvoluteFFT.cxx
2 //
3 // David Adams
4 // April 2019
5 //
6 // Test AdcDeconvoluteFFT.
7 
8 #include <string>
9 #include <iostream>
10 #include <fstream>
11 #include <sstream>
12 #include <vector>
13 #include <iomanip>
17 #include <TRandom.h>
18 #include <TH1F.h>
19 #include <TCanvas.h>
20 #include "TLegend.h"
21 
22 #undef NDEBUG
23 #include <cassert>
24 
25 using std::string;
26 using std::cout;
27 using std::endl;
28 using std::ostringstream;
29 using std::istringstream;
30 using std::ofstream;
32 using std::vector;
33 using std::setw;
34 using std::fixed;
35 
36 using Index = unsigned int;
37 using Name = std::string;
38 
39 //**********************************************************************
40 namespace {
41 
42 class Plot {
43 public:
44  TPadManipulator man;
45  TLegend* pleg;
46  Index nhst = 0;
47  Plot() : man(1000, 500) {
48  pleg = man.addLegend(0.70, 0.60, 0.90, 0.85);
49  }
50  int addHist(const AdcSignalVector& vals, Name nam, int col, int lwid, int lsty) {
51  Index nsam = vals.size();
52  Name hnam = "h" + nam;
53  Name httl = "Samples for " + nam + "; Tick; Signal";
54  TH1* ph = new TH1F(hnam.c_str(), httl.c_str(), nsam, 0, nsam);
55  ph->SetStats(0);
56  ph->SetLineWidth(lwid);
57  ph->SetLineStyle(lsty);
58  ph->SetLineColor(col);
59  for ( Index isam=0; isam<nsam; ++isam ) ph->Fill(isam+1, vals[isam]);
60  Name dopt = nhst ? "HIST SAME" : "HIST";
61  man.add(ph, dopt);
62  if ( nhst == 0 ) {
63  man.setGridY();
64  }
65  pleg->AddEntry(ph, nam.c_str(), "l");
66  ++nhst;
67  //pcan->cd();
68  //ph->Draw(dopt.c_str());
69  return 0;
70  }
71  int print(Name fname, float ymin =0.0, float ymax =0.0) {
72  if ( ymax > ymin ) man.setRangeY(ymin, ymax);
73  man.draw();
74  man.print(fname.c_str());
75  return 0;
76  }
77 };
78 
79 } // end unnamed namepsace
80 //**********************************************************************
81 
82 int test_AdcDeconvoluteFFT(bool useExistingFcl, float noiseLev, float sigmaFilter, Index nsam) {
83  const string mypre = "test_AdcDeconvoluteFFT";
84  const string myname = mypre + ": ";
85 #ifdef NDEBUG
86  cout << myname << "NDEBUG must be off." << endl;
87  abort();
88 #endif
89  string line = "-----------------------------";
90 
91  cout << myname << line << endl;
92  string fclfile = "test_AdcDeconvoluteFFT.fcl";
93  AdcSignalVector res = {0.03, 0.31, 0.77, 0.97, 1.00, 0.96, 0.92, 0.84, 0.77, 0.66,
94  0.53, 0.39, 0.31, 0.26, 0.21, 0.17, 0.14, 0.115, 0.08, 0.06,
95  0.043, 0.029, 0.022, 0.015, 0.009, 0.005, 0.003, 0.001};
96  float sum = 0.0;
97  for ( float val : res ) sum += val;
98  for ( float& val : res ) val /= sum;
99  if ( ! useExistingFcl ) {
100  cout << myname << "Creating top-level FCL." << endl;
101  ofstream fout(fclfile.c_str());
102  fout << "tools: {" << endl;
103  fout << " mydco: {" << endl;
104  fout << " tool_type: AdcDeconvoluteFFT" << endl;
105  fout << " LogLevel: 2" << endl;
106  fout << " ResponseVectors: [[";
107  for ( Index ival=0; ival<res.size(); ++ival ) {
108  if ( ival ) fout << ", ";
109  fout << res[ival];
110  }
111  fout << "]]" << endl;
112  fout << " Action: 1" << endl;
113  fout << " ResponseCenters: [5]" << endl;
114  fout << " SmoothVectors: [[-0.5, 0.5]]" << endl;
115  fout << " SmoothScales: [1.0]" << endl;
116  fout << " GausFilterSigmas: [" << sigmaFilter << "]" << endl;
117  fout << " LowFilterWidths: [-1.0]" << endl;
118  fout << " LowFilterPowers: [ 2.0]" << endl;
119  fout << " IndexMapTool: \"\"" << endl;
120  fout << " }" << endl;
121  fout << "}" << endl;
122  fout << "tools.mycondir: @local::tools.mydco" << endl;
123  fout << "tools.mycondir.Action: 4" << endl;
124  fout << "tools.myconfft: @local::tools.mydco" << endl;
125  fout << "tools.myconfft.Action: 2" << endl;
126  fout.close();
127  } else {
128  cout << myname << "Using existing top-level FCL." << endl;
129  }
130 
131  cout << myname << line << endl;
132  cout << myname << "Fetching tool manager." << endl;
134  assert ( ptm != nullptr );
135  DuneToolManager& tm = *ptm;
136  tm.print();
137  assert( tm.toolNames().size() == 3 );
138 
139  cout << myname << line << endl;
140  cout << myname << "Fetching tool." << endl;
141  auto pcondir = tm.getPrivate<TpcDataTool>("mycondir");
142  assert( pcondir != nullptr );
143  auto pconfft = tm.getPrivate<TpcDataTool>("myconfft");
144  assert( pconfft != nullptr );
145  auto pdco = tm.getPrivate<TpcDataTool>("mydco");
146  assert( pdco != nullptr );
147 
148  cout << myname << line << endl;
149  cout << myname << "Create true samples." << endl;
150  AdcSignalVector samsTru(nsam, 0.0);
151  samsTru[20] = 1000.0;
152 
153  cout << myname << line << endl;
154  cout << myname << "Create data by direct convolution with response." << endl;
155  AdcChannelData acd;
156  acd.setEventInfo(123, 456);
157  acd.setChannelInfo(789);
158  acd.samples = samsTru;
159  assert( acd.samples.size() == nsam );
160  DataMap ret = pcondir->update(acd);
161  assert( ret == 0 );
162  cout << " # samples: " << acd.samples.size() << endl;
163  cout << " # DFT amps: " << acd.dftmags.size() << endl;
164  cout << " # DFT phases: " << acd.dftphases.size() << endl;
165  assert( acd.samples.size() == nsam );
166  AdcSignalVector samsDatNoNoise = acd.samples;
167 
168  cout << myname << line << endl;
169  cout << myname << "Create data by FFT convolution with response." << endl;
170  AdcChannelData acdchk;
171  acdchk.setEventInfo(123, 456);
172  acdchk.setChannelInfo(789);
173  acdchk.samples = samsTru;
174  DataMap retchk = pconfft->update(acdchk);
175  assert( retchk == 0 );
176  cout << " # samples: " << acdchk.samples.size() << endl;
177  cout << " # DFT amps: " << acdchk.dftmags.size() << endl;
178  cout << " # DFT phases: " << acdchk.dftphases.size() << endl;
179  assert( acdchk.samples.size() == nsam );
180  AdcSignalVector samsDatCheck = acdchk.samples;
181 
182  cout << myname << line << endl;
183  cout << myname << "Compare convolutions." << endl;
184  Index nbad = 0;
185  ostringstream eout;
186  eout << " Direct FFT" << endl;
187  for ( Index isam=0; isam<nsam; ++isam ) {
188  float qdir = samsDatNoNoise[isam];
189  float qfft = samsDatCheck[isam];
190  eout << setw(5) << isam << ": " << fixed << setw(10) << qdir << fixed << setw(10) << qfft;
191  if ( fabs(qfft - qdir) > 1.e-5 ) {
192  eout << " *";
193  ++nbad;
194  }
195  eout << endl;
196  }
197  if ( nbad != 0 ) {
198  cout << myname << "Convolutions disagree:" << endl;
199  cout << eout.str() << endl;
200  }
201 
202  assert( nbad == 0 );
203 
204  cout << myname << line << endl;
205  cout << myname << "Add noise." << endl;
206  if ( noiseLev > 0.0 ) {
207  for ( float& val : acd.samples ) val += gRandom->Gaus(0.0, noiseLev);
208  }
209  AdcSignalVector samsDat = acd.samples;
210 
211  cout << myname << line << endl;
212  cout << myname << "Deconvolute." << endl;
213  ret = pdco->update(acd);
214  AdcSignalVector samsDco = acd.samples;
215  assert( ret == 0 );
216 
217  cout << myname << line << endl;
218  cout << myname << "Check integrals." << endl;
219  float ares = 0.0;
220  for ( float sig : res ) ares += sig;
221  float aTru = 0.0;
222  for ( float sig : samsTru ) aTru += sig;
223  float aDatNoNoise = 0.0;
224  for ( float sig : samsDatNoNoise ) aDatNoNoise += sig;
225  float aDat = 0.0;
226  for ( float sig : samsDat ) aDat += sig;
227  float aDco = 0.0;
228  for ( float sig : samsDco ) aDco += sig;
229  cout << myname << " Response function: " << ares << endl;
230  cout << myname << " Input signal: " << aTru << endl;
231  cout << myname << " Data without noise: " << aDatNoNoise << endl;
232  cout << myname << " Data with noise: " << aDat << endl;
233  cout << myname << " Deconvolution: " << aDco << endl;
234 
235  cout << myname << line << endl;
236  cout << myname << "Plot spectra." << endl;
237  Plot timPlot;
238  timPlot.addHist(samsTru, "True", 3, 2, 1);
239  //timPlot.addHist(res, "Response", 5, 2, 1);
240  timPlot.addHist(samsDatNoNoise, "DataNoNoise", 38, 4, 1);
241  timPlot.addHist(samsDat, "Data", 1, 2, 1);
242  timPlot.addHist(samsDco, "Deconvoluted", 2, 2, 1);
243  timPlot.print(mypre + "_time.png", -50, 200);
244 
245  cout << myname << line << endl;
246  cout << myname << "Done." << endl;
247  return 0;
248 }
249 
250 //**********************************************************************
251 
252 int main(int argc, char* argv[]) {
253  bool useExistingFcl = false;
254  float noiseLev = 4.0;
255  float sigmaFilter = 3.0;
256  Index nsam = 200;
257  if ( argc > 1 ) {
258  string sarg(argv[1]);
259  if ( sarg == "-h" ) {
260  cout << "Usage: " << argv[0] << " [KEEPFCL [NOISE [SIGMAFIL [NSAM [SETSEED]]]]]" << endl;
261  cout << " KEEPFCL = if 1 (or true), existing FCL file is used [false]" << endl;
262  cout << " NOISE = noise sigma [2.0]" << endl;
263  cout << " SIGMAFIL = Filter sigma [Tick] [2.0]" << endl;
264  cout << " NSAM = data length [200]" << endl;
265  cout << " SETSEED = if 1 (or true), a new sed is used for the noise" << endl;
266  return 0;
267  }
268  useExistingFcl = sarg == "true" || sarg == "1";
269  }
270  if ( argc > 2 ) {
271  string sarg(argv[2]);
272  istringstream ssarg(sarg);
273  ssarg >> noiseLev;
274  }
275  if ( argc > 3 ) {
276  string sarg(argv[3]);
277  istringstream ssarg(sarg);
278  ssarg >> sigmaFilter;
279  }
280  if ( argc > 4 ) {
281  string sarg(argv[4]);
282  nsam = std::stoi(sarg);
283  }
284  if ( argc > 5 ) {
285  string sarg(argv[5]);
286  if ( sarg == "true" || sarg == "1" ) gRandom->SetSeed();
287  }
288  return test_AdcDeconvoluteFFT(useExistingFcl, noiseLev, sigmaFilter, nsam);
289 }
290 
291 //**********************************************************************
int setGridY(bool flag=true)
int add(unsigned int ipad, TObject *pobj, std::string sopt="", bool replace=false)
TLegend * addLegend(double x1, double y1, double x2, double y2)
const std::vector< std::string > & toolNames() const
std::string string
Definition: nybbler.cc:12
int main(int argc, char *argv[])
struct vector vector
ChannelGroupService::Name Name
void print() const
unsigned int Index
tm
Definition: demo.py:21
void setChannelInfo(ChannelInfoPtr pchi)
const double e
void setEventInfo(EventInfoPtr pevi)
AdcSignalVector dftphases
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
std::unique_ptr< T > getPrivate(std::string name)
void line(double t, double *p, double &x, double &y, double &z)
std::vector< AdcSignal > AdcSignalVector
Definition: AdcTypes.h:22
int setRangeY(double y1, double y2)
int print(std::string fname, std::string spat="{,}")
static DuneToolManager * instance(std::string fclname="", int dbg=1)
AdcSignalVector samples
AdcSignalVector dftmags
int test_AdcDeconvoluteFFT(bool useExistingFcl, float noiseLev, float sigmaFilter, Index nsam)
QTextStream & endl(QTextStream &s)