test_fieldresp.cxx
Go to the documentation of this file.
1 #include "WireCellUtil/Testing.h"
2 
3 
4 /// needed to pretend like we are doing WCT internals
9 
10 
11 #include "TFile.h"
12 #include "TGraph.h"
13 
14 #include <iostream>
15 
16 using namespace WireCell;
17 using namespace std;
18 
19 int main(int argc, char* argv[])
20 {
21  std::string frfname = "ub-10-wnormed.json.bz2";
22  if (argc>1) {
23  frfname = argv[1];
24  cerr << "Using command line field response file: " << frfname << endl;
25  }
26  else {
27  cerr << "Using default field response file: " << frfname << endl;
28  }
29 
30 
31  /// WCT internals, normally user code does not need this
32  {
34  pm.add("WireCellSigProc");
35  auto ifrcfg = Factory::lookup<IConfigurable>("FieldResponse");
36  auto cfg = ifrcfg->default_configuration();
37  cfg["filename"] = frfname;
38  ifrcfg->configure(cfg);
39  }
40 
41  auto ifr = Factory::find<IFieldResponse>("FieldResponse");
42 
43  // Get full, "fine-grained" field responses defined at impact
44  // positions.
45  Response::Schema::FieldResponse fr = ifr->field_response();
46 
47  cerr << "FR with " << fr.planes[0].paths.size() << " responses per plane\n";
48  Assert(fr.planes[0].paths.size() == 21*6);
49 
50  // Make a new data set which is the average FR
52 
53  cerr << "FR with " << fravg.planes[0].paths.size() << " responses per plane\n";
54  /// fixme: why is this is producing 22 responses per plane and not 21?
55 
56 
57 
58  double period = 0.5 * units::microsecond;
59 
60  // convolute with electronics response function
62  WireCell::Binning tbins(Response::as_array(fravg.planes[0]).cols(), 0, Response::as_array(fravg.planes[0]).cols() * fravg.period);
64  auto ewave = ce.generate(tbins);
65  Waveform::scale(ewave, 1.2*4096/2000.);
66  elec = Waveform::dft(ewave);
67 
68  std::complex<float> fine_period(fravg.period,0);
69 
70  // TFile *file = new TFile("temp1.root","RECREATE");
71  // TGraph **gu = new TGraph*[21];
72  // TGraph **gv = new TGraph*[21];
73  // TGraph **gw = new TGraph*[21];
74  // TGraph *gtemp;
75 
76  Waveform::realseq_t wfs(9594);
77  Waveform::realseq_t ctbins(9594);
78  for (int i=0;i!=9594;i++){
79  ctbins.at(i) = i * period;
80  }
81  Waveform::realseq_t ftbins(Response::as_array(fravg.planes[0]).cols());
82  for (int i=0;i!=Response::as_array(fravg.planes[0]).cols();i++){
83  ftbins.at(i) = i * fravg.period;
84  }
85 
86  // Convert each average FR to a 2D array
87  for (int ind=0; ind<3; ++ind) {
88  auto arr = Response::as_array(fravg.planes[ind]);
89 
90  // do FFT for response ...
91  Array::array_xxc c_data = Array::dft_rc(arr,0);
92  int nrows = c_data.rows();
93  int ncols = c_data.cols();
94 
95  for (int irow = 0; irow < nrows; ++irow){
96  for (int icol = 0; icol < ncols; ++ icol){
97  c_data(irow,icol) = c_data(irow,icol) * elec.at(icol) * fine_period;
98  }
99  }
100 
101  arr = Array::idft_cr(c_data,0);
102 
103 
104  // figure out how to do fine ... shift (good ...)
105  auto arr1 = arr.block(0,0,nrows,100);
106  arr.block(0,0,nrows,ncols-100) = arr.block(0,100,nrows,ncols-100);
107  arr.block(0,ncols-100,nrows,100) = arr1;
108 
109 
110 
111  // redigitize ...
112  for (int irow = 0; irow < nrows; ++ irow){
113  // gtemp = new TGraph();
114 
115  int fcount = 1;
116  for (int i=0;i!=9594;i++){
117  double ctime = ctbins.at(i);
118 
119  if (fcount < 1000)
120  while(ctime > ftbins.at(fcount)){
121  fcount ++;
122  if (fcount >= 1000) break;
123  }
124 
125 
126  if(fcount < 1000){
127  //interpolate between fbins.at(fcount - 1) and fbins.at(fcount)
128  wfs.at(i) = (ctime - ftbins.at(fcount-1)) /fravg.period * arr(irow,fcount-1) + (ftbins.at(fcount)-ctime)/fravg.period * arr(irow,fcount);
129 
130  }else{
131  wfs.at(i) = 0;
132  }
133 
134  // gtemp->SetPoint(i,ctime/units::microsecond,wfs.at(i)/units::mV*(-1));
135  }
136 
137 
138  // if (ind ==0){
139  // gu[irow] = gtemp;
140  // // for (int icol=0; icol < ncols; ++ icol){
141  // // gu[irow]->SetPoint(icol,icol*0.1,arr(irow,icol)/units::mV*(-1));
142  // // }
143  // }else if (ind==1){
144  // gv[irow] = gtemp;
145  // // for (int icol=0; icol < ncols; ++ icol){
146  // // gv[irow]->SetPoint(icol,icol*0.1,arr(irow,icol)/units::mV*(-1));
147  // // }
148  // }else if (ind==2){
149  // gw[irow] = gtemp;
150  // // for (int icol=0; icol < ncols; ++ icol){
151  // // gw[irow]->SetPoint(icol,icol*0.1,arr(irow,icol)/units::mV*(-1));
152  // // }
153  // }
154 
155 
156  }
157 
158  //cerr << "FRavg: plane " << ind << ": " << arr.rows() << " X " << arr.cols() << " " << fravg.period/units::microsecond << endl;
159  }
160 
161  // file->cd();
162 
163  // for (int i=0;i!=21;i++){
164  // gu[i]->Write(Form("gu_%d",i));
165  // gv[i]->Write(Form("gv_%d",i));
166  // gw[i]->Write(Form("gw_%d",i));
167  // }
168  // file->Close();
169 
170 
171 
172 
173 
174 
175  // for (size_t i=0;i!=ewave.size();i++){
176  // std::cout << i *0.1 << " " << ewave.at(i) / (units::mV/units::fC) << std::endl;
177  // }
178 
179 
180 
181 }
array_xxf idft_cr(const array_xxc &arr, int dim=0)
Definition: Array.cxx:153
Sequence< real_t > realseq_t
Definition: Waveform.h:31
std::string string
Definition: nybbler.cc:12
Schema::FieldResponse wire_region_average(const Schema::FieldResponse &fr)
Definition: Response.cxx:99
Eigen::ArrayXXcf array_xxc
A complex, 2D array.
Definition: Array.h:57
STL namespace.
static const double microsecond
Definition: Units.h:94
array_xxc dft_rc(const array_xxf &arr, int dim=0)
Definition: Array.cxx:40
static const double mV
Definition: Units.h:180
cfg
Definition: dbjson.py:29
A functional object caching gain and shape.
Definition: Response.h:165
static PluginManager & instance()
Binning tbins(nticks, t0, t0+readout_time)
#define Assert
Definition: Testing.h:7
Array::array_xxf as_array(const Schema::PlaneResponse &pr)
Definition: Response.cxx:279
compseq_t dft(realseq_t seq)
Definition: Waveform.cxx:141
std::vector< PlaneResponse > planes
List of PlaneResponse objects.
Definition: Response.h:78
static const double fC
Definition: Units.h:113
Definition: Main.h:22
void scale(Sequence< Val > &seq, Val scalar)
Scale (multiply) sequence values by scalar.
Definition: Waveform.h:146
Plugin * add(const std::string &plugin_name, const std::string &libname="")
Add a plugin. If libname is not given, try to derive it.
unsigned nrows(sqlite3 *db, std::string const &tablename)
Definition: helpers.cc:84
int main(int argc, char *argv[])
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34
Hold info about multiple plane responses in the detector.
Definition: Response.h:75
QTextStream & endl(QTextStream &s)