WaveformProcessor_service.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file WaveformProcessor_service.cc
3 ///
4 /// see Calibrator/func/ADCShapeFitTable
5 /// Calibrator/art/Calibrator_service
6 ///
7 /// \author ram2aq@virginia.edu
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <iostream>
11 
13 
14 template<class T> inline T sqr(T x){return x*x;}
15 
16 namespace novaddt
17 {
18 
19  // Number of points we expect to receive, to find the best peak match with
20  const unsigned int kNumFineTimingADCPoints = 4; // Including the DCS origin
21 
22  const double kZeroOffsetSamples = -1;
23 
24  const double kSamplesPerOffset = 1./64;
25 
26  const unsigned int kNumTDCPerSample = 32;
27 
28 
29  // this method is no longer used, but could come in handy later
30  /*TDC OffsetToTDC(uint32_t const& t0, uint64_t const& high_time,
31  uint8_t const& offset)
32  {
33  TDC retval;
34  retval.val = high_time<<32 | (t0 + (uint64_t)(offset *
35  (kNumTDCPerSample * kSamplesPerOffset)));
36  retval.frac = offset%2 * 50;
37  return retval;
38  }*/
39 
40  //......................................................................
41  /// \brief Helper function for \a ADCShapeFit inner loop.
42  ///
43  /// Also used independently to calculate the best fit adc peak height (the
44  /// return value).
45  inline double GetExpectations(double t0,
46  double riseTime, double fallTime,
47  const int16_t* obs,
48  double* exps)
49  {
50  for(unsigned int t = 0; t < kNumFineTimingADCPoints; ++t){
51  // This is the ideal ASIC formula from SimpleReadout
52  exps[t] = (t > t0) ? exp(-(t-t0)/fallTime)*(1-exp(-(t-t0)/riseTime)) : 0;
53  }
54  // We can calculate the best fit baseline and normalization at this offset
55  // analytically (just write the chisq expression, and take derivatives with
56  // norm or base), like so:
57  double U, V, W, X, Z;
58  U = V = W = X = Z = 0;
59  for(unsigned int i = 0; i < kNumFineTimingADCPoints; ++i){
60  U += exps[i];
61  V += 1;
62  W += obs[i];
63  X += exps[i]*exps[i];
64  Z += obs[i]*exps[i];
65  }
66  const double norm = (W*U-V*Z)/(U*U-X*V);
67  const double base = (W*X-Z*U)/(X*V-U*U);
68 
69  // Correct all the expectations for the best guess baseline and norm
70  for(unsigned int n = 0; n < kNumFineTimingADCPoints; ++n){
71  exps[n] = norm*exps[n]+base;
72  }
73  return norm;
74  }
75 
76  /*
77 static float decodeVal8Bit(unsigned int input){
78 
79  // First get the exponent
80  // unsigned int expo = input >> 5;
81 
82  // Next get the Mantisa
83  // unsigned int mantisa = input & 0x1f;
84 
85  // The value that we return includes an averaging so that
86  // we minimize the error (but this converts it to a floating point)
87  // return (pow(2, (expo-1))*(mantisa + 32.5) - 0.5);
88 
89  // This line is the fully optimized version of this
90 
91  // Pass through values below our floor
92  //if(input < 64){return input;}
93 
94  // This computes the value if its over the pass through floor
95  // return (0x1 << ((input >> 5) -1)) *((input & 0x1f) + 32.5) - 0.5);
96 
97  // Final hand optimized version
98 
99  bool sign = input >> 8;
100  input = input & 0xff;
101 
102  if(sign)
103  return (input < 64) ? input : (0x1 << ((input >> 5) -1)) *((input & 0x1f) + 32.5) - 0.5;
104  else{
105  input = 256 - input;
106  return -((input < 64) ? input : (0x1 << ((input >> 5) -1)) *((input & 0x1f) + 32.5) - 0.5);
107  }
108 }
109 
110 */
111 
112 // Helper function that returns a mask with the Nth bit set
113 inline unsigned int bitMask(unsigned int n){
114  return (n < 32 ) ? (0x1 << n) : 0x0 ;
115 }
116 
117 
118 // Helper function that reutrns a mask with the lower N bits set
119 inline unsigned int nbitMask(unsigned int n){
120  unsigned int mask = 0xffffffff;
121  mask = mask << (32 - n);
122  mask = mask >> (32 - n);
123  return mask;
124 }
125 
126 
127 // Helper function that returns the highest bit set in an unsigned int
128 int highestBit(unsigned int input){
129  // Now determine the position of the highest non-zero bit
130  int highBit_idx=31;
131  while(highBit_idx > 0){
132  if( bitMask(highBit_idx) & input ){
133  // We have found a non-zero bit
134  // so we break out
135  return highBit_idx;
136  }
137  --highBit_idx;
138  }
139  return -1;
140 }
141 
142 
143  unsigned int NBitInt2Float(float _input, int inputBits, int outputBitsMantisa, int outputBitsExpo){
144  // First construct the bit mask for the input
145  unsigned int InputMask = 0;
146  unsigned int nBitIn = 0;
147  int highBit = 0;
148  int mantisaShift = 0;
149  unsigned int expo = 0;
150  unsigned int fp_rep=0;
151 
152  int input = (int)_input;
153  bool neg = false;
154  if(input < 0){
155  input = -1*input;
156  neg = true;
157  }
158 
159  if(input > 4095)
160  input = 4095;
161 
162  // If the intput value is smaller than the smallest exponent, just return the
163  // value without further encoding
164  if(input < (0x1 << (outputBitsMantisa+1) )){return (neg) ? 256-input: 256+input;}
165 
166  // Create an input mask to mask down the input (i.e. ensure that any random
167  // high bits above the specified length are masked down
168  InputMask = nbitMask(inputBits);
169  // p(InputMask);
170 
171  // Now mask off to the correct bits
172  nBitIn = input & InputMask;
173  // p(nBitIn);
174 
175  // Find the highest non-zero bit in the input word
176  highBit = highestBit(nBitIn);
177  if(highBit < 0){return 0;} // Check for the error code
178 
179  //printf("High Bit: %d\n",highBit);
180 
181  // Comput the shift that is required
182  mantisaShift = highBit - outputBitsMantisa;
183 
184  // printf("M Shift: %d\n",mantisaShift);
185 
186  // Shift the input word down
187  if(mantisaShift > 0){
188  nBitIn = nBitIn >> mantisaShift;
189  }
190  // p(nBitIn);
191 
192  // Now mask off the top bit
193  nBitIn &= nbitMask(outputBitsMantisa);
194 
195  // p(nBitIn);
196 
197  // Construct the exponent
198 
199  if(mantisaShift > 0){
200  expo = mantisaShift + 1;
201  expo &= nbitMask(outputBitsExpo);
202  }
203  //p(expo);
204 
205  // Shift the exponent into the correct postion
206  expo = expo << outputBitsMantisa;
207 
208  //p(expo);
209 
210  // Finally combine the mantisa and exponent
211  fp_rep = expo | nBitIn;
212 
213  //p(fp_rep);
214 
215  return (neg) ? 256 - fp_rep: 256 + fp_rep;
216 }
217 
218 
219 
220  //......................................................................
221  uint8_t ADCShapeFit(int16_t adc1, int16_t adc2, int16_t adc3,
222  double riseTime, double fallTime,
223  bool& goodTime)
224  {
225  // The observations, in chronological order
226  const int16_t obs[kNumFineTimingADCPoints] = {0, adc1, adc2, adc3};
227  double bestchisq = 1e10; // infinity
228 
229  // Sane default
230  uint8_t bestoffset = -kZeroOffsetSamples/kSamplesPerOffset;
231  goodTime = false;
232 
233  // Scan through possible offsets of the peak, looking for the best match
234  uint8_t offset = 0;
235  do{
236  // In TDC
237  const double t0 = (double(offset)*kSamplesPerOffset+kZeroOffsetSamples);
238  // Expectations
239  double exps[kNumFineTimingADCPoints];
240  const double norm = GetExpectations(t0, riseTime, fallTime, obs, exps);
241 
242  if(norm < 0) continue;
243  // Assuming equal, uncorrelated, errors on all points. (Including
244  // correlations was not found to improve resolution).
245 
246  double chisq = 0;
247  for(unsigned int i = 0; i < kNumFineTimingADCPoints; ++i)
248  chisq += sqr(obs[i]-exps[i]);
249  assert(chisq >= 0);
250 
251  if(chisq < bestchisq){
252  bestchisq = chisq;
253  bestoffset = offset;
254  goodTime = true;
255  }
256  } while(++offset != 0);
257 
258  if(bestoffset == 0 || bestoffset == 255){
259  /*
260  std::cerr << "Best time offset is at edge of range. "
261  << "If this happens repeatedly there's a "
262  << "timing problem. Details: "
263  << adc0 << " " << adc1 << " " << adc2 << std::endl;
264  */
265  bestoffset = -kZeroOffsetSamples/kSamplesPerOffset;
266  goodTime = false;
267  }
268 
269  return bestoffset;
270  }
271 
272 
273  //---------------------------------------------------------------
276  : fSampleTime (pset.get< double > ("SampleTime"))
277  , fRiseTime (pset.get< double > ("RiseTime"))
278  , fFallTime (pset.get< double > ("FallTime"))
279  , fUseMulti (pset.get< bool > ("UseMulti"))
280  , fUseMulti4ADC(pset.get< bool > ("UseMulti4ADC"))
281  , fDoHighADCFits(pset.get< bool > ("DoHighADCFits"))
282  {
283  this->reconfigure(pset);
285  }
286  //----------------------------------------------------------------
288  {
289  }
290 
291  //----------------------------------------------------------------
293  {
294  cet::search_path sp("FW_SEARCH_PATH");
295  sp.find_file(pset.get< std::string >("ShapeTableFilename"), fShapeTableFilename);
296  std::cout << pset.get< std::string > ("ShapeTableFilename") << std::endl;
297  assert(!fShapeTableFilename.empty());
298  }
299 
300  //---------------------------------------------------------------
302  {
303  if (fUseMulti && !fTableLoaded) { //load the table if it doesn't already exist
304  nTableLoads++;
305  std::cout << "Loading table from file for " << nTableLoads << " time." << std::endl;
306  TFile f(fShapeTableFilename.c_str());
307  assert(!f.IsZombie());
308 
309  TTree* tr = (TTree*)f.Get("shape_lookup");
310  assert(tr);
311 
312  const unsigned int N = tr->GetEntries();
313  assert(N == kNumTableEntries);
314 
315  uint8_t offset;
316  bool good;
317  uint16_t adc;
318  tr->SetBranchAddress("offset", &offset);
319  tr->SetBranchAddress("good", &good);
320  if (fUseMulti4ADC)
321  tr->SetBranchAddress("adc", &adc);
322 
323  for (unsigned int n=0; n<N; n++){
324  tr->GetEntry(n);
325  fTable[n] = offset;
326  fGood[n] = good;
327  if (fUseMulti4ADC)
328  fADCTable[n] = adc;
329  }
330 
331  fTableLoaded = true;
332  } // end if (!fTableLoaded)
333  }
334 
335 
336  void WaveformProcessor::process(uint32_t const& t0, uint64_t const& high_time,
337  int16_t const& adc1, int16_t const& adc2,
338  int16_t const& adc3, ADC& adcpeak,
339  bool& goodTime, TDC& retval)
340  {
341  uint8_t offset=0;
342 
343 
344  int16_t _adc1 = NBitInt2Float(adc1, 12, 5, 3);
345  int16_t _adc2 = NBitInt2Float(adc2, 12, 5, 3);
346  int16_t _adc3 = NBitInt2Float(adc3, 12, 5, 3);
347 
348  if(_adc3 >= kADC3Min && _adc3 < kADC3Max &&
349  _adc1 >= kADC1Min && _adc1 < _adc3+kADC1RelMax &&
350  _adc2 >= kADC2Min && _adc2 < kADC2Max){
351 
352  // And this is the offset it'll be at, based on aforementioned geometry
353  const int tableIdx = ((kADC2Max-kADC2Min)* // Wedge thickness
354  ((kADC3Min-kADC1Min)*(_adc3-kADC3Min+1) + // Rectangular part
355  kADC1RelMax*(_adc3-kADC3Min)+
356  ((_adc3-kADC3Min)*(_adc3-kADC3Min-1))/2+(_adc1-kADC3Min))+ // Triangular part
357  (_adc2-kADC2Min)); // Distance through wedge
358 
359  assert(tableIdx >= 0 && tableIdx < int(kNumTableEntries));
360  offset = fTable[tableIdx];
361  goodTime = fGood[tableIdx];
362 
363  }
364  else{
365  // It's not in our table, perform the fit if fDoHighADCFits is true
366  if(fDoHighADCFits)
367  offset = ADCShapeFit(adc1, adc2, adc3, fRiseTime, fFallTime, goodTime);
368  else{
369  offset = 0;
370  goodTime = false;
371  }
372  }
373 
374  // We need to extract the best normalization, so repeat the relevant part
375  // of the inner loop at that offset.
376  // --- For speed the normalization has been saved in a different table,
377  // --- the use of which is set as a fhicl parameter.
378  if(goodTime){
379  //Perform the normalization:
380  // The observations, in chronological order
381  //const int16_t obs[kNumFineTimingADCPoints] = {0, adc1, adc2, adc3};
382  //const double t0 = (double(offset)*kSamplesPerOffset+kZeroOffsetSamples);
383  //double junk[kNumFineTimingADCPoints]; // Expectations
384  //const double norm = GetExpectations(t0, fRiseTime, fFallTime, obs, junk);
385 
386  // The value the peak of the curve would have before scaling by norm
387  // Analytic calculation from the functional form in ReadoutSim
388  //static const double peak = pow(fRiseTime/(fRiseTime+fFallTime),
389  // fRiseTime/fFallTime)-pow(fRiseTime/(fRiseTime+fFallTime),
390  // 1+fRiseTime/fFallTime);
391  //adcpeak = norm*peak;
392 
393  if (fUseMulti4ADC) { //Extract from lookup table for speed:
394  // And this is the offset it'll be at, based on aforementioned geometry
395  const int tableIdx = ((kADC2Max-kADC2Min)* // Wedge thickness
396  ((kADC3Min-kADC1Min)*(_adc3-kADC3Min+1) + // Rectangular part
397  kADC1RelMax*(_adc3-kADC3Min)+
398  ((_adc3-kADC3Min)*(_adc3-kADC3Min-1))/2+(_adc1-kADC3Min))+ // Triangular part
399  (_adc2-kADC2Min)); // Distance through wedge
400 
401 
402 
403  adcpeak = fADCTable[tableIdx];
404  } // else case handled in initmulti()
405 
406  }
407  else{
408  // Unless the fit failed, in which case we're probably safest just doing
409  // this.
410  adcpeak = ADC(adc3);
411  }
412 
413  retval.val = (high_time<<32 | (t0 + (uint64_t)(offset *
414  (kNumTDCPerSample * kSamplesPerOffset))));
415  retval.frac = (uint8_t)(offset%2 * 50);
416  }
417 
418  void WaveformProcessor::initmulti() { // initialize TDC and ADC members
419  //using multipoint
420  //check that multipoint can be used
421  const uint32_t version = fNano->getHeader()->getVersion();
422 
424 
425  const unsigned int nSamples = conv.getNSamples(version);
426  const unsigned int nPretrig = conv.getNPretriggeredSamples(version);
427 
428  if ( (nSamples < 3) || // too few samples
429  (nSamples <= nPretrig) || //missing peak
430  (version == 0) ) //non-multipoint
431  { initdcs(); return; }
432 
433 
434  TDC tempTDC;
435  ADC tempADC;
436  uint32_t t0 = fNano->getTimeStamp();
437  uint64_t high_time = (uint64_t) fMicro->getHighWord();
438  int16_t adc1 = (int16_t)fNano->getValue(1) - (int16_t)fNano->getValue(0);
439  int16_t adc2 = (int16_t)fNano->getValue(2) - (int16_t)fNano->getValue(0);
440  int16_t adc3 = (int16_t)fNano->getValue(3) - (int16_t)fNano->getValue(0);
441  bool goodTime = true;
442  process(t0, high_time, adc1, adc2, adc3, tempADC, goodTime, tempTDC);
443 
444  if(adc3 < 0) adc3 = 0; // Check if adc3 is negative, ADC is unsigned
445  if(tempADC.val < 0) tempADC = (ADC)0; // Check if adc3 is negative, ADC is unsigned
446  if(!fUseMulti4ADC) tempADC = (ADC)adc3;
447 
448  fTDC = tempTDC;
449  fADC = tempADC;
450  }
451 
452  void WaveformProcessor::setup(uint32_t const& t0, uint64_t const& high_time,
453  int16_t const& adc1, int16_t const& adc2, int16_t const& adc3,
454  int16_t const& dcs)
455  {
456  // For use by RawDigit2DAQHit in the offline.
457  // dcs is identical to adc3 for 4-point readout
458  if (!fUseMulti) {
459  fTDC.val = high_time<<32 | t0;
460  fADC = ADC(dcs);
461  return;
462  }
463  bool goodTime = true;
464  TDC tempTDC;
465  ADC tempADC;
466  process(t0, high_time, adc1, adc2, adc3, tempADC, goodTime, tempTDC);
467  if (!goodTime) { // something's gone wrong with the fit, use DCS
468  fTDC.val = high_time<<32 | t0;
469  fADC = ADC(adc3);
470  return;
471  }
472  else {
473  fTDC = tempTDC;
474  if (fUseMulti4ADC)
475  fADC = tempADC;
476  else
477  fADC = ADC(adc3);
478 
479  }
480  }
481 
482  void WaveformProcessor::initdcs() { // initialize TDC and ADC members
483  //using dcs
484 
485  //auto trash = fMicro->getFloatingMicroSlice();
486  fTDC.val = ((uint64_t)fMicro->getHighWord() << 32)
487  | fNano->getTimeStamp();
488 
489  fTDC.frac = 0;
490  int16_t adc3 = (int16_t)fNano->getValue(3) - (int16_t)fNano->getValue(0);
491  if(adc3 > 0)
492  fADC = (ADC)adc3;
493  else
494  fADC = (ADC)0;
495 
496  }
497 
499 
500 } // end namespace novaddt
WaveformProcessor(fhicl::ParameterSet const &pset, art::ActivityRegistry &reg)
static const int kADC1Min
IMPLEMENT_MAIN_STANDARD IMPLEMENT_MAIN_setBufferSource daqdataformats::RawNanoSliceHeader * getHeader() const
Get the NanoSlice header pointer.
Definition: RawNanoSlice.h:57
static const int kADC2Min
fraction_type frac
Definition: BaseProducts.h:35
value_type val
Definition: BaseProducts.h:34
static const int kADC3Min
#define DEFINE_ART_SERVICE(svc)
Definition: ServiceMacros.h:93
Float_t x1[n_points_granero]
Definition: compare.C:5
const double kSamplesPerOffset
void process(uint32_t const &t0, uint64_t const &high_time, int16_t const &adc1, int16_t const &adc2, int16_t const &adc3, ADC &adcpeak, bool &goodTime, TDC &retval)
const unsigned int kNumFineTimingADCPoints
daqdataformats::RawNanoSlice * fNano
unsigned int bitMask(unsigned int n)
int highestBit(unsigned int input)
GlobalSignal< detail::SignalResponseType::LIFO, void(Run const &)> sPostBeginRun
std::bitset< kNumTableEntries > fGood
void reconfigure(const fhicl::ParameterSet &pset)
std::string find_file(std::string const &filename) const
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
static const int kADC1RelMax
Definition: Run.h:31
Float_t Z
Definition: plot.C:38
double GetExpectations(double t0, double riseTime, double fallTime, const int16_t *obs, double *exps)
Helper function for ADCShapeFit inner loop.
uint32_t getNSamples(const version_t ver) const
Get the number of samples.
uint32_t getNPretriggeredSamples(const version_t ver) const
Get number of pretriggered samples.
static const unsigned int kNumTableEntries
T get(std::string const &key) const
Definition: ParameterSet.h:231
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const double kZeroOffsetSamples
std::array< uint16_t, kNumTableEntries > fADCTable
static const int kADC2Max
std::array< uint8_t, kNumTableEntries > fTable
Definition: run.py:1
void setup(daqdataformats::RawMicroSlice *theMicro, daqdataformats::RawNanoSlice *theNano)
OStream cout
Definition: OStream.cxx:6
const XML_Char * version
Definition: expat.h:187
daqdataformats::RawMicroSlice * fMicro
Service for extracting timing and pulse height information from traces with multipoint readout...
unsigned int nTableLoads
unsigned int NBitInt2Float(float _input, int inputBits, int outputBitsMantisa, int outputBitsExpo)
Float_t norm
unsigned int nbitMask(unsigned int n)
assert(nhit_max >=nhit_nbins)
value_type val
Definition: BaseProducts.h:65
void postBeginRun(art::Run const &run)
const unsigned int kNumTDCPerSample
static const int kADC3Max
double T
Definition: Xdiff_gwt.C:5
uint32_t getHighWord() const
Get method for lower half of timing word.
Float_t X
Definition: plot.C:38
#define W(x)
uint8_t ADCShapeFit(int16_t adc1, int16_t adc2, int16_t adc3, double riseTime, double fallTime, bool &goodTime)
enum BeamMode string