Public Member Functions | Private Attributes | List of all members
WireCell::Gen::ImpactTransform Class Reference

#include <ImpactTransform.h>

Public Member Functions

 ImpactTransform (IPlaneImpactResponse::pointer pir, BinnedDiffusion_transform &bd)
 
virtual ~ImpactTransform ()
 
Waveform::realseq_t waveform (int wire) const
 

Private Attributes

IPlaneImpactResponse::pointer m_pir
 
BinnedDiffusion_transform & m_bd
 
int m_num_group
 
int m_num_pad_wire
 
std::vector< std::map< int, IImpactResponse::pointer > > m_vec_map_resp
 
std::vector< std::vector< std::tuple< int, int, double > > > m_vec_vec_charge
 
std::vector< int > m_vec_impact
 
Array::array_xxf m_decon_data
 
int m_start_ch
 
int m_end_ch
 
int m_start_tick
 
int m_end_tick
 
Log::logptr_t log
 

Detailed Description

An ImpactTransform transforms charge on impact positions into waveforms via 2D FFT.

Definition at line 18 of file ImpactTransform.h.

Constructor & Destructor Documentation

Gen::ImpactTransform::ImpactTransform ( IPlaneImpactResponse::pointer  pir,
BinnedDiffusion_transform &  bd 
)

Definition at line 9 of file ImpactTransform.cxx.

10  :m_pir(pir), m_bd(bd)
11  , log(Log::logger("sim"))
12 {
13 
14  // for (int i=0;i!=210;i++){
15  // double pos = -31.5 + 0.3*i+1e-9;0
16  // m_pir->closest(pos);
17  // }
18 
19  // arrange the field response (210 in total, pitch_range/impact)
20  // number of wires nwires ...
21  m_num_group = std::round(m_pir->pitch()/m_pir->impact())+1; // 11
22  m_num_pad_wire = std::round((m_pir->nwires()-1)/2.); // 10
23 
24  const auto pimpos = m_bd.pimpos();
25  // const int nsamples = m_bd.tbins().nbins();
26  //const auto rb = pimpos.region_binning();
27  //const int nwires = rb.nbins();
28 
29  //
30 
31  //std::cout << m_num_group << " " << m_num_pad_wire << std::endl;
32  for (int i=0;i!=m_num_group;i++){
33  double rel_cen_imp_pos ;
34  if (i!=m_num_group-1){
35  rel_cen_imp_pos = -m_pir->pitch()/2.+m_pir->impact()*i+1e-9;
36  }else{
37  rel_cen_imp_pos = -m_pir->pitch()/2.+m_pir->impact()*i-1e-9;
38  }
39  m_vec_impact.push_back(std::round(rel_cen_imp_pos/m_pir->impact()));
40  std::map<int, IImpactResponse::pointer> map_resp; // already in freq domain
41 
42  for (int j=0;j!=m_pir->nwires();j++){
43  map_resp[j-m_num_pad_wire] = m_pir->closest(rel_cen_imp_pos - (j-m_num_pad_wire)*m_pir->pitch());
44  Waveform::compseq_t response_spectrum = map_resp[j-m_num_pad_wire]->spectrum();
45 
46  // std::cout << i << " " << j << " " << rel_cen_imp_pos - (j-m_num_pad_wire)*m_pir->pitch()<< " " << response_spectrum.size() << std::endl;
47  }
48  //std::cout << m_vec_impact.back() << std::endl;
49  // std::cout << rel_cen_imp_pos << std::endl;
50  // std::cout << map_resp.size() << std::endl;
51  m_vec_map_resp.push_back(map_resp);
52 
53  //Eigen::SparseMatrix<float> *mat = new Eigen::SparseMatrix<float>(nsamples,nwires);
54  // mat.reserve(Eigen::VectorXi::Constant(nwires,1000));
55  //m_vec_spmatrix.push_back(mat);
56 
57  std::vector<std::tuple<int,int, double> > vec_charge; // ch, time, charge
58  m_vec_vec_charge.push_back(vec_charge);
59  }
60 
61  // m_bd.get_charge_matrix(m_vec_spmatrix, m_vec_impact);
62  //std::cout << nwires << " " << nsamples << std::endl;
63 
64 
65 
66  // now work on the charge part ...
67  // trying to sampling ...
68  m_bd.get_charge_vec(m_vec_vec_charge, m_vec_impact);
69  //std::cout << nwires << " " << nsamples << std::endl;
70 
71  // for (size_t i=0;i!=m_vec_vec_charge.size();i++){
72  // std::cout << m_vec_vec_charge[i].size() << std::endl;
73  // }
74 
75  // length and width ...
76 
77 
78  //
79 
80 
81  // std::cout << nwires << " " << nsamples << std::endl;
82  std::pair<int,int> impact_range = m_bd.impact_bin_range(m_bd.get_nsigma());
83  std::pair<int,int> time_range = m_bd.time_bin_range(m_bd.get_nsigma());
84 
85  // std::cout << impact_range.first << " " << impact_range.second << " " << time_range.first << " " << time_range.second << std::endl;
86 
87  int start_ch = std::floor(impact_range.first*1.0/(m_num_group-1))-1;
88  int end_ch = std::ceil(impact_range.second*1.0/(m_num_group-1))+2;
89  if ( (end_ch-start_ch)%2==1) end_ch += 1;
90  int start_tick = time_range.first-1;
91  int end_tick = time_range.second+2;
92  if ( (end_tick-start_tick)%2==1 ) end_tick += 1;
93 
94 
95  // m_decon_data = Array::array_xxf::Zero(nwires+2*m_num_pad_wire,nsamples);
96  // for saving the accumulated wire data in the time frequency domain ...
97  // adding no padding now, it make the FFT slower, need some other methods ...
98 
99  int npad_wire =0;
100  const size_t ntotal_wires = fft_best_length(end_ch - start_ch + 2 * m_num_pad_wire,1);
101 
102  // pow(2,std::ceil(log(end_ch - start_ch + 2 * m_num_pad_wire)/log(2)));
103  // if (nwires == 2400){
104  // if (ntotal_wires > 2500)
105  // ntotal_wires = 2500;
106  // }else if (nwires ==3456){
107  // if (ntotal_wires > 3600)
108  // ntotal_wires = 3600;
109  // npad_wire=72; //3600
110  //}
111  npad_wire= (ntotal_wires -end_ch + start_ch)/2;
112  m_start_ch = start_ch - npad_wire;
113  m_end_ch = end_ch + npad_wire;
114  //std::cout << start_ch << " " << end_ch << " " << npad_wire << " " << start_tick << " " << end_tick << " " << m_start_ch << " " << m_end_ch << std::endl;
115 
116 
117  int npad_time = m_pir->closest(0)->waveform_pad();
118  const size_t ntotal_ticks = fft_best_length(end_tick - start_tick + npad_time);
119 
120  // pow(2,std::ceil(log(end_tick - start_tick + npad_time)/log(2)));
121  // if (ntotal_ticks >9800 && nsamples <9800 && nsamples >9550)
122  // ntotal_ticks = 9800;
123  npad_time = ntotal_ticks - end_tick + start_tick;
124  m_start_tick = start_tick;
125  m_end_tick = end_tick + npad_time;
126 
127  // m_end_tick = 16384;//nsamples;
128  // m_start_tick = 0;
129  // // std::cout << m_start_tick << " " << m_end_tick << std::endl;
130  // int npad_time = 0;
131  // int ntotal_ticks = pow(2,std::ceil(log(nsamples + npad_time)/log(2)));
132  // if (ntotal_ticks >9800 && nsamples <9800)
133  // ntotal_ticks = 9800
134  // npad_time = ntotal_ticks - nsamples;
135  // m_start_tick = 0;
136  // m_end_tick = ntotal_ticks;
137 
138 
139  Array::array_xxc acc_data_f_w = Array::array_xxc::Zero(end_ch-start_ch+2*npad_wire, m_end_tick - m_start_tick);
140 
141  int num_double = (m_vec_vec_charge.size()-1)/2;
142  //int num_double = (m_vec_spmatrix.size()-1)/2;
143 
144  // speed up version , first five
145  for (int i=0;i!=num_double;i++){
146  //if (i!=0) continue;
147  // std::cout << i << std::endl;
148  Array::array_xxc c_data = Array::array_xxc::Zero(end_ch-start_ch+2*npad_wire,m_end_tick - m_start_tick);
149 
150  // fill normal order
151  for (size_t j=0;j!=m_vec_vec_charge.at(i).size();j++){
152  c_data(std::get<0>(m_vec_vec_charge.at(i).at(j))+npad_wire-start_ch,std::get<1>(m_vec_vec_charge.at(i).at(j))-m_start_tick) += std::get<2>(m_vec_vec_charge.at(i).at(j));
153  }
154  // std::cout << i << " " << m_vec_vec_charge.at(i).size() << std::endl;
155  m_vec_vec_charge.at(i).clear();
156  m_vec_vec_charge.at(i).shrink_to_fit();
157 
158  // useing matrix form ...
159  // for (int k=0; k<m_vec_spmatrix.at(i)->outerSize(); ++k)
160  // for (Eigen::SparseMatrix<float>::InnerIterator it(*m_vec_spmatrix.at(i),k); it; ++it){
161  // c_data(it.col()+npad_wire-start_ch,it.row()-m_start_tick) = it.value();
162  // }
163  // delete m_vec_spmatrix.at(i);
164  // //m_vec_spmatrix.at(i).setZero();
165  // //m_vec_spmatrix.at(i).resize(0,0);
166 
167  // fill reverse order
168  int ii=num_double*2-i;
169  for (size_t j=0;j!=m_vec_vec_charge.at(ii).size();j++){
170  c_data(end_ch+npad_wire-1-std::get<0>(m_vec_vec_charge.at(ii).at(j)),std::get<1>(m_vec_vec_charge.at(ii).at(j))-m_start_tick) += std::complex<float>(0,std::get<2>(m_vec_vec_charge.at(ii).at(j)));
171  }
172  // std::cout << ii << " " << m_vec_vec_charge.at(ii).size() << std::endl;
173  m_vec_vec_charge.at(ii).clear();
174  m_vec_vec_charge.at(ii).shrink_to_fit();
175  // for (int k=0; k<m_vec_spmatrix.at(ii)->outerSize(); ++k)
176  // for (Eigen::SparseMatrix<float>::InnerIterator it(*m_vec_spmatrix.at(ii),k); it; ++it){
177  // c_data(it.col()+npad_wire-start_ch,it.row()-m_start_tick) = it.value();
178  // }
179  // delete m_vec_spmatrix.at(ii);
180  // // m_vec_spmatrix.at(ii).setZero();
181  // //m_vec_spmatrix.at(ii).resize(0,0);
182 
183 
184  // Do FFT on time
185  c_data = Array::dft_cc(c_data,0);
186  // Do FFT on wire
187  c_data = Array::dft_cc(c_data,1);
188 
189  // std::cout << i << std::endl;
190  {
191  Array::array_xxc resp_f_w = Array::array_xxc::Zero(end_ch-start_ch+2*npad_wire,m_end_tick - m_start_tick);
192  {
193  Waveform::compseq_t rs1 = m_vec_map_resp.at(i)[0]->spectrum();
194  // do a inverse FFT
195  Waveform::realseq_t rs1_t = Waveform::idft(rs1);
196  // pick the first xxx ticks
197  Waveform::realseq_t rs1_reduced(m_end_tick-m_start_tick,0);
198  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
199  if (icol >= int(rs1_t.size())) break;
200  rs1_reduced.at(icol) = rs1_t[icol];
201  }
202  // do a FFT
203  rs1 = Waveform::dft(rs1_reduced);
204 
205  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
206  resp_f_w(0,icol) = rs1[icol];
207  }
208  }
209 
210 
211  for (int irow = 0; irow!=m_num_pad_wire;irow++){
212  Waveform::compseq_t rs1 = m_vec_map_resp.at(i)[irow+1]->spectrum();
213  Waveform::realseq_t rs1_t = Waveform::idft(rs1);
214  Waveform::realseq_t rs1_reduced(m_end_tick-m_start_tick,0);
215  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
216  if (icol >= int(rs1_t.size())) break;
217  rs1_reduced.at(icol) = rs1_t[icol];
218  }
219  rs1 = Waveform::dft(rs1_reduced);
220  Waveform::compseq_t rs2 = m_vec_map_resp.at(i)[-irow-1]->spectrum();
221  Waveform::realseq_t rs2_t = Waveform::idft(rs2);
222  Waveform::realseq_t rs2_reduced(m_end_tick-m_start_tick,0);
223  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
224  if (icol >= int(rs2_t.size())) break;
225  rs2_reduced.at(icol) = rs2_t[icol];
226  }
227  rs2 = Waveform::dft(rs2_reduced);
228  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
229  resp_f_w(irow+1,icol) = rs1[icol];
230  resp_f_w(end_ch-start_ch-1-irow+2*npad_wire,icol) = rs2[icol];
231  }
232  }
233  //std::cout << i << std::endl;
234 
235  // Do FFT on wire for response // slight larger
236  resp_f_w = Array::dft_cc(resp_f_w,1); // Now becomes the f and f in both time and wire domain ...
237  // multiply them together
238  c_data = c_data * resp_f_w;
239  }
240 
241 
242  // Do inverse FFT on wire
243  c_data = Array::idft_cc(c_data,1);
244 
245  // Add to wire result in frequency
246  acc_data_f_w += c_data;
247 
248  }
249 
250  // std::cout << "ABC : " << std::endl;
251 
252  // central region ...
253  {
254  int i = num_double;
255  // fill response array in frequency domain
256 
257  Array::array_xxc data_f_w;
258  {
259  Array::array_xxf data_t_w = Array::array_xxf::Zero(end_ch-start_ch+2*npad_wire,m_end_tick-m_start_tick);
260  // fill charge array in time-wire domain // slightly larger
261  for (size_t j=0;j!=m_vec_vec_charge.at(i).size();j++){
262  data_t_w(std::get<0>(m_vec_vec_charge.at(i).at(j))+npad_wire-start_ch,std::get<1>(m_vec_vec_charge.at(i).at(j))-m_start_tick) += std::get<2>(m_vec_vec_charge.at(i).at(j));
263  // std::cout << std::get<1>(m_vec_vec_charge.at(i).at(j)) << std::endl;
264  }
265  // std::cout << i << " " << m_vec_vec_charge.at(i).size() << std::endl;
266  m_vec_vec_charge.at(i).clear();
267  m_vec_vec_charge.at(i).shrink_to_fit();
268  // for (int k=0; k<m_vec_spmatrix.at(i)->outerSize(); ++k)
269  // for (Eigen::SparseMatrix<float>::InnerIterator it(*m_vec_spmatrix.at(i),k); it; ++it){
270  // data_t_w(it.col()+npad_wire-start_ch,it.row()-m_start_tick) = it.value();
271  // }
272  // delete m_vec_spmatrix.at(i);
273  // // m_vec_spmatrix.at(i).setZero();
274  // // m_vec_spmatrix.at(i).resize(0,0);
275 
276 
277 
278  // Do FFT on time
279  data_f_w = Array::dft_rc(data_t_w,0);
280  // Do FFT on wire
281  data_f_w = Array::dft_cc(data_f_w,1);
282  }
283 
284  {
285  Array::array_xxc resp_f_w = Array::array_xxc::Zero(end_ch-start_ch+2*npad_wire,m_end_tick-m_start_tick);
286 
287  {
288  Waveform::compseq_t rs1 = m_vec_map_resp.at(i)[0]->spectrum();
289  // Array::array_xxc temp_resp_f_w = Array::array_xxc::Zero(2*m_num_pad_wire+1,nsamples);
290  // for (int icol = 0; icol != nsamples; icol++){
291  // temp_resp_f_w(0,icol) = rs1[icol];
292  // }
293  // Array::array_xxf temp_resp_t_w = Array::idft_cr(temp_resp_f_w,0).block(0,0,2*m_num_pad_wire+1,m_end_tick-m_start_tick);
294  // temp_resp_f_w = Array::dft_rc(temp_resp_t_w,0);
295 
296  // do a inverse FFT
297  Waveform::realseq_t rs1_t = Waveform::idft(rs1);
298  // pick the first xxx ticks
299  Waveform::realseq_t rs1_reduced(m_end_tick-m_start_tick,0);
300  // std::cout << rs1.size() << " " << nsamples << " " << m_end_tick << " " << m_start_tick << std::endl;
301  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
302  if (icol >= int(rs1_t.size())) break;
303  rs1_reduced.at(icol) = rs1_t[icol];
304  // std::cout << icol << " " << rs1_t[icol] << std::endl;
305  }
306  // do a FFT
307  rs1 = Waveform::dft(rs1_reduced);
308 
309  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
310  // std::cout << icol << " " << rs1[icol] << " " << temp_resp_f_w(0,icol) << std::endl;
311  resp_f_w(0,icol) = rs1[icol];
312  }
313  }
314  for (int irow = 0; irow!=m_num_pad_wire;irow++){
315  Waveform::compseq_t rs1 = m_vec_map_resp.at(i)[irow+1]->spectrum();
316  Waveform::realseq_t rs1_t = Waveform::idft(rs1);
317  Waveform::realseq_t rs1_reduced(m_end_tick-m_start_tick,0);
318  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
319  if (icol >= int(rs1_t.size())) break;
320  rs1_reduced.at(icol) = rs1_t[icol];
321  }
322  rs1 = Waveform::dft(rs1_reduced);
323  Waveform::compseq_t rs2 = m_vec_map_resp.at(i)[-irow-1]->spectrum();
324  Waveform::realseq_t rs2_t = Waveform::idft(rs2);
325  Waveform::realseq_t rs2_reduced(m_end_tick-m_start_tick,0);
326  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
327  if (icol >= int(rs2_t.size())) break;
328  rs2_reduced.at(icol) = rs2_t[icol];
329  }
330  rs2 = Waveform::dft(rs2_reduced);
331  for (int icol = 0; icol != m_end_tick-m_start_tick; icol++){
332  resp_f_w(irow+1,icol) = rs1[icol];
333  resp_f_w(end_ch-start_ch-1-irow+2*npad_wire,icol) = rs2[icol];
334  }
335  // for (int icol = 0; icol != nsamples; icol++){
336  // resp_f_w(irow+1,icol) = rs1[icol];
337  // resp_f_w(end_ch-start_ch-1-irow+2*npad_wire,icol) = rs2[icol];
338  // }
339  }
340  // Do FFT on wire for response // slight larger
341  resp_f_w = Array::dft_cc(resp_f_w,1); // Now becomes the f and f in both time and wire domain ...
342  // multiply them together
343  data_f_w = data_f_w * resp_f_w;
344  }
345 
346  // Do inverse FFT on wire
347  data_f_w = Array::idft_cc(data_f_w,1);
348 
349  // Add to wire result in frequency
350  acc_data_f_w += data_f_w;
351  }
352 
353  //m_decon_data = Array::array_xxc::Zero(nwires,nsamples);
354  // if (npad_wire!=0){
355  acc_data_f_w = Array::idft_cc(acc_data_f_w,0);//.block(npad_wire,0,nwires,nsamples);
356  Array::array_xxf real_m_decon_data = acc_data_f_w.real();
357  Array::array_xxf img_m_decon_data = acc_data_f_w.imag().colwise().reverse();
358  m_decon_data = real_m_decon_data + img_m_decon_data;
359 
360  // std::cout << real_m_decon_data(40,5182) << " " << img_m_decon_data(40,5182) << std::endl;
361  // std::cout << real_m_decon_data(40,5182-m_start_tick) << " " << img_m_decon_data(40,5182-m_start_tick) << std::endl;
362 
363  //}else{
364  // Array::array_xxc temp_m_decon_data = Array::idft_cc(acc_data_f_w,0);
365  // Array::array_xxf real_m_decon_data = temp_m_decon_data.real();
366  // Array::array_xxf img_m_decon_data = temp_m_decon_data.imag().rowwise().reverse();
367  // m_decon_data = real_m_decon_data + img_m_decon_data;
368  // }
369 
370  // // prepare FFT, loop 11 of them ... (older version)
371  // for (size_t i=0;i!=m_vec_vec_charge.size();i++){
372  // // fill response array in frequency domain
373  // if (i!=10) continue;
374 
375  // Array::array_xxc data_f_w;
376  // {
377  // Array::array_xxf data_t_w = Array::array_xxf::Zero(nwires+2*npad_wire,nsamples);
378  // // fill charge array in time-wire domain // slightly larger
379  // for (size_t j=0;j!=m_vec_vec_charge.at(i).size();j++){
380  // data_t_w(std::get<0>(m_vec_vec_charge.at(i).at(j))+npad_wire,std::get<1>(m_vec_vec_charge.at(i).at(j))) += std::get<2>(m_vec_vec_charge.at(i).at(j));
381  // }
382  // m_vec_vec_charge.at(i).clear();
383 
384  // // Do FFT on time
385  // data_f_w = Array::dft_rc(data_t_w,0);
386  // // Do FFT on wire
387  // data_f_w = Array::dft_cc(data_f_w,1);
388  // }
389 
390  // {
391  // Array::array_xxc resp_f_w = Array::array_xxc::Zero(nwires+2*npad_wire,nsamples);
392  // {
393  // Waveform::compseq_t rs1 = m_vec_map_resp.at(i)[0]->spectrum();
394  // for (int icol = 0; icol != nsamples; icol++){
395  // resp_f_w(0,icol) = rs1[icol];
396  // }
397  // }
398  // for (int irow = 0; irow!=m_num_pad_wire;irow++){
399  // Waveform::compseq_t rs1 = m_vec_map_resp.at(i)[irow+1]->spectrum();
400  // Waveform::compseq_t rs2 = m_vec_map_resp.at(i)[-irow-1]->spectrum();
401  // for (int icol = 0; icol != nsamples; icol++){
402  // resp_f_w(irow+1,icol) = rs1[icol];
403  // resp_f_w(nwires-1-irow+2*npad_wire,icol) = rs2[icol];
404  // }
405  // }
406  // // Do FFT on wire for response // slight larger
407  // resp_f_w = Array::dft_cc(resp_f_w,1); // Now becomes the f and f in both time and wire domain ...
408  // // multiply them together
409  // data_f_w = data_f_w * resp_f_w;
410  // }
411 
412  // // Do inverse FFT on wire
413  // data_f_w = Array::idft_cc(data_f_w,1);
414 
415  // // Add to wire result in frequency
416  // acc_data_f_w += data_f_w;
417  // }
418  // m_vec_vec_charge.clear();
419 
420  // // do inverse FFT on time for the final results ...
421 
422  // if (npad_wire!=0){
423  // Array::array_xxf temp_m_decon_data = Array::idft_cr(acc_data_f_w,0);
424  // m_decon_data = temp_m_decon_data.block(npad_wire,0,nwires,nsamples);
425  // }else{
426  // m_decon_data = Array::idft_cr(acc_data_f_w,0);
427  // }
428 
429 
430  // std::cout << m_decon_data(40,5195-m_start_tick)/units::mV << " " << m_decon_data(40,5195-m_start_tick)/units::mV << std::endl;
431 
432  // m_vec_spmatrix.clear();
433  //m_vec_spmatrix.shrink_to_fit();
434 
435  // int nrows = resp_f_w.rows();
436  // int ncols = resp_f_w.cols();
437  log->debug("ImpactTransform: # of channels: {} # of ticks: {}",
438  m_decon_data.rows(), m_decon_data.cols());
439 
440 } // constructor
WireCell::Gen::ImpactTransform::m_pir
IPlaneImpactResponse::pointer m_pir
Definition: ImpactTransform.h:20
WireCell::Array::dft_cc
array_xxc dft_cc(const array_xxc &arr, int dim=1)
Definition: Array.cxx:67
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
WireCell::Gen::BinnedDiffusion_transform::time_bin_range
std::pair< int, int > time_bin_range(double nsigma=0.0) const
Definition: BinnedDiffusion_transform.cxx:456
WireCell::fft_best_length
std::size_t fft_best_length(size_t nsamples, bool keep_odd_even=false)
WireCell::Array::array_xxc
Eigen::ArrayXXcf array_xxc
A complex, 2D array.
Definition: Array.h:57
WireCell::Array::dft_rc
array_xxc dft_rc(const array_xxf &arr, int dim=0)
Definition: Array.cxx:40
WireCell::Gen::BinnedDiffusion_transform::get_nsigma
double get_nsigma() const
Definition: BinnedDiffusion_transform.h:101
WireCell::Waveform::dft
compseq_t dft(realseq_t seq)
Definition: Waveform.cxx:141
e
const double e
Definition: gUpMuFluxGen.cxx:165
WireCell::Array::idft_cc
array_xxc idft_cc(const array_xxc &arr, int dim=1)
Definition: Array.cxx:125
WireCell::Gen::ImpactTransform::m_vec_impact
std::vector< int > m_vec_impact
Definition: ImpactTransform.h:30
WireCell::Gen::ImpactTransform::m_bd
BinnedDiffusion_transform & m_bd
Definition: ImpactTransform.h:21
WireCell::Log::logger
logptr_t logger(std::string name)
Definition: Logging.cxx:71
WireCell::Gen::BinnedDiffusion_transform::pimpos
const Pimpos & pimpos() const
Definition: BinnedDiffusion_transform.h:57
WireCell::Waveform::idft
realseq_t idft(compseq_t spec)
Definition: Waveform.cxx:149
pimpos
Pimpos pimpos(nwires, min_wire_pitch, max_wire_pitch)
WireCell::Gen::ImpactTransform::m_vec_vec_charge
std::vector< std::vector< std::tuple< int, int, double > > > m_vec_vec_charge
Definition: ImpactTransform.h:26
WireCell::Gen::ImpactTransform::m_vec_map_resp
std::vector< std::map< int, IImpactResponse::pointer > > m_vec_map_resp
Definition: ImpactTransform.h:25
WireCell::Gen::BinnedDiffusion_transform::get_charge_vec
void get_charge_vec(std::vector< std::vector< std::tuple< int, int, double > > > &vec_vec_charge, std::vector< int > &vec_impact)
Definition: BinnedDiffusion_transform.cxx:223
WireCell::Gen::BinnedDiffusion_transform::impact_bin_range
std::pair< int, int > impact_bin_range(double nsigma=0.0) const
Definition: BinnedDiffusion_transform.cxx:439
WireCell::Waveform::compseq_t
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34
WireCell::Array::array_xxf
Eigen::ArrayXXf array_xxf
A real, 2D array.
Definition: Array.h:54
Gen::ImpactTransform::~ImpactTransform ( )
virtual

Definition at line 444 of file ImpactTransform.cxx.

445 {
446 
447 }

Member Function Documentation

Waveform::realseq_t Gen::ImpactTransform::waveform ( int  wire) const

Return the wire's waveform. If the response functions are just field response (ie, instantaneous current) then the waveforms are expressed as current integrated over each sample bin and thus in units of charge. If the response functions include electronics response then the waveforms are in units of voltage representing the sampling of the output of the FEE amplifiers.

Definition at line 450 of file ImpactTransform.cxx.

451 {
452  const int nsamples = m_bd.tbins().nbins();
453  if (iwire < m_start_ch || iwire >= m_end_ch){
454  return Waveform::realseq_t(nsamples, 0.0);
455  }else{
456  Waveform::realseq_t wf(nsamples, 0.0);
457  for (int i=0;i!=nsamples;i++){
458  if (i>=m_start_tick && i < m_end_tick){
459  wf.at(i) = m_decon_data(iwire-m_start_ch,i-m_start_tick);
460  }else{
461  //wf.at(i) = 1e-25;
462  }
463  //std::cout << m_decon_data(iwire-m_start_ch,i-m_start_tick) << std::endl;
464  }
465 
466  if (m_pir->closest(0)->long_aux_waveform().size()>0){
467  // now convolute with the long-range response ...
468  const size_t nlength = fft_best_length(nsamples + m_pir->closest(0)->long_aux_waveform_pad());
469 
470  //nlength = nsamples;
471 
472  // std::cout << nlength << " " << nsamples + m_pir->closest(0)->long_aux_waveform_pad() << std::endl;
473 
474  wf.resize(nlength,0);
475  Waveform::realseq_t long_resp = m_pir->closest(0)->long_aux_waveform();
476  long_resp.resize(nlength,0);
477  Waveform::compseq_t spec = Waveform::dft(wf);
478  Waveform::compseq_t long_spec = Waveform::dft(long_resp);
479  for (size_t i=0;i!=nlength;i++){
480  spec.at(i) *= long_spec.at(i);
481  }
482  wf = Waveform::idft(spec);
483  wf.resize(nsamples,0);
484  }
485 
486  return wf;
487  }
488 
489 }
WireCell::Gen::ImpactTransform::m_pir
IPlaneImpactResponse::pointer m_pir
Definition: ImpactTransform.h:20
WireCell::Waveform::realseq_t
Sequence< real_t > realseq_t
Definition: Waveform.h:31
WireCell::fft_best_length
std::size_t fft_best_length(size_t nsamples, bool keep_odd_even=false)
WireCell::Waveform::dft
compseq_t dft(realseq_t seq)
Definition: Waveform.cxx:141
WireCell::Gen::ImpactTransform::m_bd
BinnedDiffusion_transform & m_bd
Definition: ImpactTransform.h:21
WireCell::Waveform::idft
realseq_t idft(compseq_t spec)
Definition: Waveform.cxx:149
WireCell::Gen::BinnedDiffusion_transform::tbins
const Binning & tbins() const
Definition: BinnedDiffusion_transform.h:58
WireCell::Binning::nbins
int nbins() const
Definition: Binning.h:42
WireCell::Waveform::compseq_t
Sequence< complex_t > compseq_t
A complex-valued sequence, eg for discrete spectrum powers.
Definition: Waveform.h:34

Member Data Documentation

Log::logptr_t WireCell::Gen::ImpactTransform::log
private

Definition at line 37 of file ImpactTransform.h.

BinnedDiffusion_transform& WireCell::Gen::ImpactTransform::m_bd
private

Definition at line 21 of file ImpactTransform.h.

Array::array_xxf WireCell::Gen::ImpactTransform::m_decon_data
private

Definition at line 31 of file ImpactTransform.h.

int WireCell::Gen::ImpactTransform::m_end_ch
private

Definition at line 33 of file ImpactTransform.h.

int WireCell::Gen::ImpactTransform::m_end_tick
private

Definition at line 35 of file ImpactTransform.h.

int WireCell::Gen::ImpactTransform::m_num_group
private

Definition at line 23 of file ImpactTransform.h.

int WireCell::Gen::ImpactTransform::m_num_pad_wire
private

Definition at line 24 of file ImpactTransform.h.

IPlaneImpactResponse::pointer WireCell::Gen::ImpactTransform::m_pir
private

Definition at line 20 of file ImpactTransform.h.

int WireCell::Gen::ImpactTransform::m_start_ch
private

Definition at line 32 of file ImpactTransform.h.

int WireCell::Gen::ImpactTransform::m_start_tick
private

Definition at line 34 of file ImpactTransform.h.

std::vector<int> WireCell::Gen::ImpactTransform::m_vec_impact
private

Definition at line 30 of file ImpactTransform.h.

std::vector<std::map<int, IImpactResponse::pointer> > WireCell::Gen::ImpactTransform::m_vec_map_resp
private

Definition at line 25 of file ImpactTransform.h.

std::vector<std::vector<std::tuple<int,int,double> > > WireCell::Gen::ImpactTransform::m_vec_vec_charge
private

Definition at line 26 of file ImpactTransform.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11