Public Types | Public Member Functions | Private Types | Private Attributes | List of all members
Tpc2dDeconvolute Class Reference

#include <Tpc2dDeconvolute.h>

Inheritance diagram for Tpc2dDeconvolute:
TpcDataTool

Public Types

using Index = unsigned int
 

Public Member Functions

 Tpc2dDeconvolute (fhicl::ParameterSet const &ps)
 
 ~Tpc2dDeconvolute () override=default
 
DataMap updateTpcData (TpcData &) const override
 
DataMap viewTpcData (const TpcData &) const override
 
Fw2dFFT & fft () const
 

Private Types

using Name = std::string
 
using IndexVector = std::vector< Index >
 
using ResponseVector = AdcSignalVector
 
using ResponseVectorVector = std::vector< ResponseVector >
 
- Private Types inherited from AdcChannelTool
using Index = unsigned int
 

Private Attributes

int m_LogLevel
 
ResponseVectorVector m_ResponseVectors
 
int m_ResponseCenter
 
Index m_FftSize
 
Name m_InPath
 
Name m_OutPath
 
float m_SampleSigma
 
float m_ChannelSigma
 
float m_LowFilterWidth
 
float m_LowFilterPower
 
std::unique_ptr< Fw2dFFT > m_pfft
 

Additional Inherited Members

- Private Member Functions inherited from TpcDataTool
virtual int forwardTpcData () const
 
- Private Member Functions inherited from AdcChannelTool
virtual ~AdcChannelTool ()=default
 
virtual DataMap update (AdcChannelData &) const
 
virtual DataMap view (const AdcChannelData &acd) const
 
virtual DataMap updateMap (AdcChannelDataMap &acds) const
 
virtual DataMap viewMap (const AdcChannelDataMap &acds) const
 
virtual bool updateWithView () const
 
virtual bool viewWithUpdate () const
 
virtual DataMap beginEvent (const DuneEventInfo &) const
 
virtual DataMap endEvent (const DuneEventInfo &) const
 
virtual DataMap close (const DataMap *dmin=nullptr)
 
- Static Private Member Functions inherited from AdcChannelTool
static int interfaceNotImplemented ()
 

Detailed Description

Definition at line 72 of file Tpc2dDeconvolute.h.

Member Typedef Documentation

using Tpc2dDeconvolute::Index = unsigned int

Definition at line 76 of file Tpc2dDeconvolute.h.

using Tpc2dDeconvolute::IndexVector = std::vector<Index>
private

Definition at line 94 of file Tpc2dDeconvolute.h.

using Tpc2dDeconvolute::Name = std::string
private

Definition at line 93 of file Tpc2dDeconvolute.h.

using Tpc2dDeconvolute::ResponseVector = AdcSignalVector
private

Definition at line 95 of file Tpc2dDeconvolute.h.

using Tpc2dDeconvolute::ResponseVectorVector = std::vector<ResponseVector>
private

Definition at line 96 of file Tpc2dDeconvolute.h.

Constructor & Destructor Documentation

Tpc2dDeconvolute::Tpc2dDeconvolute ( fhicl::ParameterSet const &  ps)

Definition at line 47 of file Tpc2dDeconvolute_tool.cc.

48 : m_LogLevel(ps.get<int>("LogLevel")),
49  m_ResponseVectors(ps.get<ResponseVectorVector>("ResponseVectors")),
50  m_ResponseCenter(ps.get<int>("ResponseCenter")),
51  m_FftSize(ps.get<Index>("FftSize")),
52  m_InPath(ps.get<Name>("InPath")),
53  m_OutPath(ps.get<Name>("OutPath")),
54  m_SampleSigma(ps.get<float>("SampleSigma")),
55  m_ChannelSigma(ps.get<float>("ChannelSigma")),
56  m_LowFilterWidth(ps.get<float>("LowFilterWidth")),
57  m_LowFilterPower(ps.get<float>("LowFilterPower")),
58  m_pfft(new Fw2dFFT(m_FftSize, 1)) {
59  const string myname = "Tpc2dDeconvolute::ctor: ";
60  if ( m_LogLevel >= 1 ) {
61  std::ios cout_state(nullptr);
62  cout_state.copyfmt(std::cout);
63  cout << myname << "Parameters:" << endl;
64  cout << myname << " LogLevel: " << m_LogLevel << endl;
65  cout << myname << " ResponseVectors: [";
66  Index indent = 22;
67  bool first = true;
68  for ( const ResponseVector& vec : m_ResponseVectors ) {
69  if ( first ) first = false;
70  else cout << ", ";
71  for ( Index isam=0; isam<vec.size(); ++isam ) {
72  float val = vec[isam];
73  if ( isam != 0 ) cout << ", ";
74  if ( isam == 0 ) cout << "\n" << myname << setw(indent) << "[";
75  else if ( 10*(isam/10) == isam ) cout << "\n" << myname << setw(indent) << " ";
76  cout << setw(11) << fixed << setprecision(7) << val;
77  }
78  cout << "]";
79  }
80  cout << endl;
81  cout << myname << setw(indent-2) << "]" << endl;
82  cout << myname << " ResponseCenter: " << m_ResponseCenter << endl;
83  cout << myname << " FftSize: " << m_FftSize << endl;
84  cout << myname << " InPath: " << m_InPath << endl;
85  cout << myname << " OutPath: " << m_OutPath << endl;
86  cout << myname << " SampleSigma: " << m_SampleSigma << endl;
87  cout << myname << " ChannelSigma: " << m_ChannelSigma << endl;
88  cout << myname << " LowFilterPower: " << m_LowFilterPower << endl;
89  cout << myname << " LowFilterWidth: " << m_LowFilterWidth << endl;
90  }
91 }
Name
ChannelGroupService::Name Name
Definition: FixedChannelGroupService_service.cc:19
Index
unsigned int Index
Definition: sspmapmaker_v1.c:14
setprecision
Q_EXPORT QTSManip setprecision(int p)
Definition: qtextstream.h:343
Tpc2dDeconvolute::ResponseVectorVector
std::vector< ResponseVector > ResponseVectorVector
Definition: Tpc2dDeconvolute.h:96
Tpc2dDeconvolute::m_ResponseVectors
ResponseVectorVector m_ResponseVectors
Definition: Tpc2dDeconvolute.h:100
genie::units::ps
static constexpr double ps
Definition: Units.h:99
test_nxjs.indent
indent
Definition: test_nxjs.py:18
setw
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
Tpc2dDeconvolute::ResponseVector
AdcSignalVector ResponseVector
Definition: Tpc2dDeconvolute.h:95
Tpc2dDeconvolute::m_pfft
std::unique_ptr< Fw2dFFT > m_pfft
Definition: Tpc2dDeconvolute.h:111
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
Tpc2dDeconvolute::~Tpc2dDeconvolute ( )
overridedefault

Member Function Documentation

Fw2dFFT & Tpc2dDeconvolute::fft ( ) const

Definition at line 293 of file Tpc2dDeconvolute_tool.cc.

293  {
294  return *m_pfft;
295 }
Tpc2dDeconvolute::m_pfft
std::unique_ptr< Fw2dFFT > m_pfft
Definition: Tpc2dDeconvolute.h:111
DataMap Tpc2dDeconvolute::updateTpcData ( TpcData &  tpd) const
overridevirtual

Reimplemented from TpcDataTool.

Definition at line 95 of file Tpc2dDeconvolute_tool.cc.

95  {
96  const string myname = "Tpc2dDeconvolute::updateTpcData: ";
97  DataMap ret;
98  // Find input data.
99  TpcData* ptpdIn = tpd.getTpcData(m_InPath);
100  if ( ptpdIn == nullptr ) {
101  cout << myname << "WARNING: Input TPC data not found at " << m_InPath << "." << endl;
102  return ret.setStatus(1);
103  }
104  TpcData::Tpc2dRoiVector& inrois = ptpdIn->get2dRois();
105  // Create output data object.
106  TpcData::Tpc2dRoiVector* poutrois = nullptr;
107  if ( m_OutPath.size() && m_OutPath != "." ) {
108  TpcData* ptpdOut = tpd.getTpcData(m_OutPath);
109  if ( ptpdOut == nullptr ) {
110  if ( m_LogLevel >= 2 ) {
111  cout << myname << "Adding data path " << m_OutPath << "." << endl;
112  }
113  ptpdOut = tpd.addTpcData(m_OutPath);
114  if ( ptpdOut == nullptr ) {
115  cout << myname << "ERROR: Unable to create output TPC data at " << m_OutPath << "." << endl;
116  return ret.setStatus(2);
117  }
118  }
119  poutrois = &ptpdOut->get2dRois();
120  if ( poutrois->size() ) {
121  // Might want to add config flag to disable this message.
122  cout << myname << "WARNING: Output 2D ROI container is not empty." << endl;
123  }
124  }
125  // Define DFT normalization.
126  RealDftNormalization dftNorm(12);
127  // Loop over 2d ROIs.
128  Index xflog = 0;
129  Index nroi = 0;
130  for ( Tpc2dRoi& roi : inrois ) {
131  // Check fft size.
132  Index ncha = roi.channelSize();
133  Index nsam = roi.sampleSize();
134  IndexArray ndxs = {ncha, nsam};
135  if ( fft().checkDataSize(ndxs) ) {
136  cout << myname << "ERROR: Skipping too-large 2D ROI: " << ncha << " x " << nsam << endl;
137  continue;
138  }
139  // Get ROI DFT.
140  if ( roi.dft() == nullptr ) {
141  const RealData& rdat = roi.data();
142  DftData* pdft = new DftData(dftNorm, ndxs);
143  fft().fftForward(rdat, *pdft, xflog);
144  roi.resetDft(pdft);
145  }
146  if ( roi.dft() == nullptr ) {
147  cout << myname << "ERROR: Forward FFT of ROI failed." << endl;
148  continue;
149  }
150  DftData& roiDft = *roi.dft();
151  // Build ncha x nsam response matrix.
152  RealData resDat(ndxs);
153  ResponseVector empty;
154  // ... Response is periodic. Choose tick offset in first period.
155  // ... Response is symmetric in channel.
156  int nsamInt = nsam;
157  int samoffInt = m_ResponseCenter;
158  while ( samoffInt < 0 ) samoffInt += nsamInt;
159  while ( samoffInt >= nsamInt ) samoffInt -= nsamInt;
160  Index samoff = samoffInt;
161  for ( Index dcha1=0; dcha1<(ncha+1)/2; ++dcha1 ) {
162  const ResponseVector& res1d =
163  dcha1 < m_ResponseVectors.size() ? m_ResponseVectors[dcha1] : empty;
164  IndexVector dchas(1, dcha1);
165  if ( dcha1 ) {
166  Index dcha2 = ncha - dcha1;
167  if ( dcha2 != dcha1 ) dchas.push_back(dcha2);
168  }
169  for ( Index dsam=0; dsam<nsam; ++dsam ) {
170  Index isam = dsam + samoff;
171  if ( isam >= nsam ) isam -= nsam;
172  float rval = isam < res1d.size() ? res1d[isam] : 0.0;
173  for ( Index dcha : dchas ) {
174  IndexArray idxs = {dcha, dsam};
175  Index idat = resDat.setValue(idxs, rval);
176  if ( m_LogLevel >= 3 ) {
177  cout << myname << " res[" << dcha << ", " << dsam << "]: " << rval << endl;
178  }
179  // This shouldn't happen...
180  if ( idat >= resDat.size() ) {
181  cout << myname << "ERROR: Unexpected invalid response indices." << endl;
182  }
183  }
184  }
185  }
186  // Get response DFT.
187  RealDftNormalization respNorm(RealDftNormalization::convolutionNormalization());
188  DftData resDft(respNorm, ndxs);
189  fft().fftForward(resDat, resDft, xflog);
190  // Build sample filter.
191  // User supplies Gaussian sigma in sample (tick). Xform is a Gaussian
192  // with sigma_freq = Nsam/(2 pi sigma_time).
193  std::vector<double> samFilt(nsam, 1.0);
194  if ( m_SampleSigma > 0.0 ) {
195  float freqSigma = nsam/(2.0*acos(-1.0)*m_SampleSigma);
196  for ( Index ifrq=0; ifrq<nsam; ++ifrq ) {
197  float xf = ifrq/freqSigma;
198  float fac = exp(-xf*xf/2.0);
199  samFilt[ifrq] = fac;
200  if ( m_LogLevel >= 3 ) cout << myname << " samfil[" << ifrq << "]: " << fac << endl;
201  }
202  }
203  // Build channel filter.
204  // User supplies Gaussian sigma in channel. Xform is a Gaussian
205  // with sigma_freq = Ncha/(2 pi sigma_time).
206  std::vector<double> chaFilt(ncha, 1.0);
207  if ( m_ChannelSigma > 0.0 ) {
208  float freqSigma = nsam/(2.0*acos(-1.0)*m_ChannelSigma);
209  for ( Index ifrq=0; ifrq<ncha; ++ifrq ) {
210  float xf = ifrq/freqSigma;
211  float fac = exp(-xf*xf/2.0);
212  chaFilt[ifrq] = fac;
213  if ( m_LogLevel >= 3 ) cout << myname << " chafil[" << ifrq << "]: " << fac << endl;
214  }
215  }
216  std::vector<double> samFiltLow(nsam, 1.0);
217  if ( m_LowFilterPower > 0 ) {
218  if ( m_LowFilterWidth == 0.0 ) {
219  samFiltLow[0] = 0.0;
220  } else {
221  samFiltLow[0] = 0.0;
222  double fac = nsam/m_LowFilterWidth;
223  for ( Index ifrq=1; ifrq<nsam; ++ifrq ) {
224  double arg = pow(fac/ifrq, m_LowFilterPower);
225  samFiltLow[ifrq] = 1.0/(1.0 + arg);
226  }
227  }
228  }
229  // Evaluate deconvoluted DFT.
230  // Copy the input ROI DFT, multiply by filter and divide by response.
231  DftData* poutDft = new DftData(roiDft);
232  bool logEntry = m_LogLevel>=4;
233  for ( Index idat=0; idat<poutDft->size(); ++idat ) {
234  auto idxs = poutDft->indexArrays(idat);
235  Index dcha = idxs[0][0];
236  Index dsam = idxs[0][1];
237  Tpc2dRoi::Dft::Complex& val = poutDft->data()[idat];
238  // For the filters, use the smaller of the wave number and its conjugate.
239  Index dchaFilt = std::min(dcha, ncha-dcha);
240  Index dsamFilt = std::min(dsam, nsam-dsam);
241  if ( logEntry ) {
242  cout << myname << "Decon[" << dcha << "," << dsam << "]";
243  cout << "[" << dchaFilt << "," << dsamFilt << "]";
244  cout << ": " << comstr(val);
245  }
246  val *= chaFilt[dchaFilt];
247  val *= samFilt[dsamFilt];
248  val *= samFiltLow[dsamFilt];
249  Tpc2dRoi::Dft::Complex resp = resDft.data()[idat];
250  if ( resp != 0.0 ) val /= resp;
251  if ( logEntry ) {
252  cout << " * " << chaFilt[dchaFilt];
253  cout << " * " << samFilt[dsamFilt];
254  cout << " / " << comstr(resp) << " = " << comstr(val) << endl;
255  }
256  }
257  if ( m_LogLevel >=3 ) {
258  cout << myname << " Input power: " << roiDft.power() << endl;
259  cout << myname << "Response power: " << resDft.power() << endl;
260  cout << myname << " Output power: " << poutDft->power() << endl;
261  }
262  // Create the output ROI and attach the DFT data.
263  Tpc2dRoi* poutroi = &roi;
264  if ( poutrois != nullptr ) {
265  poutrois->emplace_back(ncha, nsam, roi.channelOffset(), roi.sampleOffset());
266  poutroi = &poutrois->back();
267  }
268  poutroi->resetDft(poutDft);
269  // Evaluate deconvoluted reponse.
270  int dstat = fft().fftBackward(*poutroi->dft(), poutroi->data());
271  if ( dstat ) {
272  cout << myname << "ERROR: Reverse xform failed with error " << dstat << endl;
273  continue;
274  }
275  ++nroi;
276  }
277  ret.setInt("dcoNroiIn", inrois.size());
278  ret.setInt("dcoNroiOut", nroi);
279  return ret;
280 }
Real2dData::data
const std::vector< F > & data() const
Definition: Real2dData.h:117
IndexVector
std::vector< Index > IndexVector
Definition: test_FclIndexVectorMapTool.cxx:26
Fw2dFFT::fftForward
int fftForward(const Data &dat, DFT &dft, Index logLevel=0)
Definition: Fw2dFFT.cxx:93
DataMap::setStatus
DataMap & setStatus(int stat)
Definition: DataMap.h:130
FftwReal2dDftData::indexArrays
IndexArrayVector indexArrays(Index idat) const
Definition: FftwReal2dDftData.h:132
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
TpcData::addTpcData
TpcData * addTpcData(Name nam, bool copyAdcData=true)
Definition: TpcData.cxx:24
IndexArray
DftData::IndexArray IndexArray
Definition: test_Fw2dFFT.cxx:35
Index
unsigned int Index
Definition: sspmapmaker_v1.c:14
RealDftNormalization::convolutionNormalization
static RealDftNormalization convolutionNormalization()
Definition: RealDftNormalization.h:162
Fw2dFFT::fftBackward
int fftBackward(const DFT &dft, Data &dat, Index logLevel=0)
Definition: Fw2dFFT.cxx:123
Tpc2dDeconvolute::fft
Fw2dFFT & fft() const
Definition: Tpc2dDeconvolute_tool.cc:293
DftData
Tpc2dRoi::Dft DftData
Definition: Tpc2dDeconvolute_tool.cc:24
Tpc2dRoi::resetDft
void resetDft(Dft *pdft)
Definition: Tpc2dRoi.h:73
Tpc2dDeconvolute::m_ResponseVectors
ResponseVectorVector m_ResponseVectors
Definition: Tpc2dDeconvolute.h:100
DataMap::setInt
void setInt(Name name, int val)
Definition: DataMap.h:131
FftwReal2dDftData::size
Index size() const
Definition: FftwReal2dDftData.h:93
Real2dDftData::Complex
std::complex< Float > Complex
Definition: Real2dDftData.h:29
TpcData::getTpcData
TpcData * getTpcData(Name nam)
Definition: TpcData.cxx:42
cet::sqlite::min
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
TpcData::Tpc2dRoiVector
std::vector< Tpc2dRoi > Tpc2dRoiVector
Definition: TpcData.h:31
Tpc2dDeconvolute::ResponseVector
AdcSignalVector ResponseVector
Definition: Tpc2dDeconvolute.h:95
Tpc2dRoi::data
const DataArray & data() const
Definition: Tpc2dRoi.h:50
Real2dDftData::power
virtual F power(const IndexArray &ifrqs) const
Definition: Real2dDftData.h:91
TpcData::get2dRois
Tpc2dRoiVector & get2dRois()
Definition: TpcData.h:57
util::empty
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:97
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
Tpc2dRoi::dft
Dft * dft()
Definition: Tpc2dRoi.h:68
DataMap Tpc2dDeconvolute::viewTpcData ( const TpcData &  ) const
overridevirtual

Reimplemented from TpcDataTool.

Definition at line 284 of file Tpc2dDeconvolute_tool.cc.

284  {
285  const string myname = "Tpc2dDeconvolute::updateTpcData: ";
286  cout << myname << "ERROR: View of TPC data is not supported here." << endl;
287  DataMap ret;
288  return ret.setStatus(1);
289 }
DataMap::setStatus
DataMap & setStatus(int stat)
Definition: DataMap.h:130
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030

Member Data Documentation

float Tpc2dDeconvolute::m_ChannelSigma
private

Definition at line 106 of file Tpc2dDeconvolute.h.

Index Tpc2dDeconvolute::m_FftSize
private

Definition at line 102 of file Tpc2dDeconvolute.h.

Name Tpc2dDeconvolute::m_InPath
private

Definition at line 103 of file Tpc2dDeconvolute.h.

int Tpc2dDeconvolute::m_LogLevel
private

Definition at line 99 of file Tpc2dDeconvolute.h.

float Tpc2dDeconvolute::m_LowFilterPower
private

Definition at line 108 of file Tpc2dDeconvolute.h.

float Tpc2dDeconvolute::m_LowFilterWidth
private

Definition at line 107 of file Tpc2dDeconvolute.h.

Name Tpc2dDeconvolute::m_OutPath
private

Definition at line 104 of file Tpc2dDeconvolute.h.

std::unique_ptr<Fw2dFFT> Tpc2dDeconvolute::m_pfft
private

Definition at line 111 of file Tpc2dDeconvolute.h.

int Tpc2dDeconvolute::m_ResponseCenter
private

Definition at line 101 of file Tpc2dDeconvolute.h.

ResponseVectorVector Tpc2dDeconvolute::m_ResponseVectors
private

Definition at line 100 of file Tpc2dDeconvolute.h.

float Tpc2dDeconvolute::m_SampleSigma
private

Definition at line 105 of file Tpc2dDeconvolute.h.

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