ADCShapeFit.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Function to perform, and table to cache, timing fits to ADC values
3 /// \author bckhouse@caltech.edu
4 /// \date
5 ////////////////////////////////////////////////////////////////////////
6 
7 // Please, do not introduce any art dependencies into this file. It needs to be
8 // able to be compiled standalone (as part of genADCShapeFitLookupTable) to
9 // produce the lookup table file in the first place. It should be possible to
10 // do any required artyness in Calibrator.cxx
11 
13 
14 #include <cassert>
15 #include <cmath>
16 #include <iostream>
17 
18 #include "lz4hc.h"
19 
20 #include "TFile.h"
21 #include "TTree.h"
22 
23 template<class T> inline T sqr(T x){return x*x;}
24 
25 namespace calib
26 {
27  // Size of the chunks into which the compressed table is split. This is a
28  // tradeoff between speed and memory. A larger number (larger chunks) gives
29  // better compression, but slower retrieval from the table. A smaller number
30  // is faster but less memory efficient.
31  //
32  // Chunk size | Memory Runtime (FD cosmics maxopt)
33  // -------------+--------------------
34  // Uncompressed | 93MB 3.0sec/evt
35  // 1,000 | 25MB 3.1sec/evt
36  // 2,000 | 20MB 3.2sec/evt
37  // 4,000 | 17MB 3.4sec/evt
38  // 10,000 | 16MB 4.3sec/evt
39  const int kChunk = 2000;
40 
41  static_assert(ADCShapeFitTable::kNumTableEntries % kChunk == 0,
42  "Compression chunks must divide table evenly");
43 
44  // Number of points we expect to receive, to find the best peak match with
45  const unsigned int kNumFineTimingADCPoints = 4; // Including the DCS origin
46 
47  // Empirically, the best fit is never found outside this range. +3 is an
48  // upper limit on the positive side because it means the pulse starts beyond
49  // the last ADC observation.
50  const double kZeroOffsetSamples = -1;
51  // Have the capacity to consider 65,536 offsets, which is way too many to
52  // give good performance.
53  const int kOffsetStep = 16;
54  const double kSamplesPerOffset = 1./(256*64);
55  // With 4096 steps that means the maximum offset is +3 and step size is
56  // 0.49ns (in the FD).
57 
58  //......................................................................
59  /// \brief Helper function for \a ADCShapeFit inner loop.
60  ///
61  /// Also used independently to calculate the best fit adc peak height (the
62  /// return value).
63 
64  inline double GetExpectations(double t0,
65  double riseTime, double fallTime, double preAmp,
66  const int16_t* obs,
67  double* exps, double& base, int mode)
68  {
69  // There are only a few thousand distinct inputs to this function per
70  // combination of riseTime, fallTime and preAmp. Cache all possible
71  // outputs. This saves about 40% of the run time of CalHit for ND and FD
72  // jobs, for which this cache needs to be built only once. For TB jobs,
73  // two kinds of FEBs are present, so the cache is rebuilt many times as we
74  // see hits from each kind. This is *still* faster than not having a
75  // cache, although obviously a little dumb since we could save both sets of
76  // cached data with a little more work.
77 
78  // The calculation of cache_size is slightly fragile because
79  // kSamplesPerOffset and kZeroOffsetSamples are defined as doubles even
80  // though they are representing integer concepts (also the meaning of
81  // kSamplesPerOffset is the reciprocal of its name). But I think the rest
82  // of the code also breaks if kZeroOffsetSamples isn't an integer or if
83  // kSamplesPerOffset isn't the reciprocal of an integer.
84  const int steps_per_offset = 1 / kSamplesPerOffset / kOffsetStep;
85  const int num_offsets = (kNumFineTimingADCPoints - 1) - kZeroOffsetSamples;
86  const int cache_size = steps_per_offset * num_offsets;
87 
88  static double cache_ideal[cache_size] = {0};
89  static double cache_3param[cache_size] = {0};
90 
91  static double cache_riseTime = -1, cache_fallTime = -1, cache_preAmp = -1;
92 
93  if(cache_ideal[0] == 0 || riseTime != cache_riseTime ||
94  fallTime != cache_fallTime || preAmp != cache_preAmp){
95  for(int i = 0; i < cache_size; i++){
96  const double t_minus_t0 = double(i+1)/steps_per_offset;
97 
98  // This is the ideal ASIC formula from SimpleReadout
99  cache_ideal[i] = exp(-t_minus_t0/fallTime)*(1-exp(-t_minus_t0/riseTime));
100 
101  //new three parameters function
102  cache_3param[i] = (fallTime * preAmp)*(
103  (exp(-t_minus_t0/fallTime) / ((preAmp-fallTime)*(fallTime-riseTime))) -
104  (exp(-t_minus_t0/preAmp) / ((preAmp-fallTime)*(preAmp - riseTime))) +
105  (exp(-t_minus_t0/riseTime) / ((preAmp-riseTime)*(riseTime-fallTime))) );
106  }
107  }
108 
109  cache_riseTime = riseTime;
110  cache_fallTime = fallTime;
111  cache_preAmp = preAmp;
112 
113  for(unsigned int t = 0; t < kNumFineTimingADCPoints; ++t){
114  if(t <= t0){
115  exps[t] = 0;
116  }
117  else{
118  const int cachei = (t-t0) * steps_per_offset - 1;
119  exps[t] = mode == 0?cache_ideal[cachei]: cache_3param[cachei];
120  }
121  }
122 
123  // We can calculate the best fit baseline and normalization at this offset
124  // analytically (just write the chisq expression, and take derivatives with
125  // norm or base), like so:
126  double U, V, W, X, Z;
127  U = V = W = X = Z = 0;
128  for(unsigned int i = 0; i < kNumFineTimingADCPoints; ++i){
129  U += exps[i];
130  V += 1;
131  W += obs[i];
132  X += exps[i]*exps[i];
133  Z += obs[i]*exps[i];
134  }
135  const double norm = (W*U-V*Z)/(U*U-X*V);
136  base = (W*X-Z*U)/(X*V-U*U);
137 
138  // Correct all the expectations for the best guess baseline and norm
139  for(unsigned int n = 0; n < kNumFineTimingADCPoints; ++n){
140  exps[n] = norm*exps[n]+base;
141  }
142 
143  return norm;
144  }
145 
146  //......................................................................
147  uint16_t ADCShapeFit(int16_t adc1, int16_t adc2, int16_t adc3,
148  double riseTime, double fallTime, double preAmp,
149  bool& goodTime, int fMode)
150  {
151  // The observations, in chronological order
152  const int16_t obs[kNumFineTimingADCPoints] = {0, adc1, adc2, adc3};
153 
154  double bestchisq = 1e10; // infinity
155  // Sane default
156  uint16_t bestoffset = -kZeroOffsetSamples/kSamplesPerOffset;
157  goodTime = false;
158 
159  // Scan through possible offsets of the peak, looking for the best match
160  uint16_t offset = 0;
161  do{
162  // In TDC
163  const double t0 = (double(offset)*kSamplesPerOffset+kZeroOffsetSamples);
164 
165  // Expectations
166  double exps[kNumFineTimingADCPoints];
167  double junk;
168  const double norm = GetExpectations(t0, riseTime, fallTime, preAmp, obs, exps, junk, fMode);
169 
170  if(norm < 0) continue;
171 
172  // Assuming equal, uncorrelated, errors on all points. (Including
173  // correlations was not found to improve resolution).
174  double chisq = 0;
175  for(unsigned int i = 0; i < kNumFineTimingADCPoints; ++i)
176  chisq += sqr(obs[i]-exps[i]);
177  assert(chisq >= 0);
178 
179  if(chisq < bestchisq){
180  bestchisq = chisq;
181  bestoffset = offset;
182  goodTime = true;
183  }
184  } while((offset += kOffsetStep) != 0);
185 
186  if(bestoffset == 0 || bestoffset == 256*256-kOffsetStep){
187  /*
188  std::cerr << "Best time offset is at edge of range. "
189  << "If this happens repeatedly there's a "
190  << "timing problem. Details: "
191  << adc0 << " " << adc1 << " " << adc2 << std::endl;
192  */
193  bestoffset = -kZeroOffsetSamples/kSamplesPerOffset;
194  goodTime = false;
195  }
196 
197  return bestoffset;
198  }
199 
200 #ifndef ADCSHAPEFIT_FUNC_ONLY
201  //......................................................................
202  ADCShapeFitTable::ADCShapeFitTable(const std::string& fname, bool isND, bool isMC, int mode)
203  : fIsND(isND), fIsMC(isMC), fMode(mode), fNHit(0), fNMiss(0)
204  {
205  TFile f(fname.c_str());
206  assert(!f.IsZombie());
207 
208  TTree* tr = (TTree*)f.Get("shape_lookup");
209  assert(tr);
210 
211  const unsigned int N = tr->GetEntries();
212 
213  assert(N == kNumTableEntries);
214 
215  // Offsets are all two bytes to save space, which is far more resolution
216  // ADCShapeFit needs.
217  uint16_t offset;
218  bool good;
219  tr->SetBranchAddress("offset", &offset);
220  tr->SetBranchAddress("good", &good);
221 
222  std::vector<uint16_t> table(N);
223 
224  std::cout << "Loading ADCShapeFit table from " << fname;
225  for(unsigned int n = 0; n < N; ++n){
226  if(n%(N/8) == 0) std::cout << "." << std::flush;
227 
228  tr->GetEntry(n);
229  table[n] = offset;
230 
231  if(!good) fIsBad.insert(n);
232  }
233  std::cout << std::endl;
234 
235  // The maximum length a compressed block can be
236  const unsigned long bufferLen = LZ4_compressBound(kChunk*sizeof(uint16_t));
237  std::vector<char> buffer(bufferLen);
238 
239  std::cout << "Compressing table";
240  int compressedSize = 0;
241  for(unsigned int i = 0; i < kNumTableEntries/kChunk; ++i){
242  if(i%(kNumTableEntries/kChunk/8) == 0) std::cout << "." << std::flush;
243 
244  // Compress one kChunk-sized chunk of the table into buffer
245  const int outputLen = LZ4_compress_HC((char*)&table[i*kChunk],
246  &buffer.front(),
247  kChunk*sizeof(uint16_t),
248  bufferLen,
250 
251  assert(outputLen > 0);
252 
253  // Store the result into the i'th entry of fCompresssedTable
254  fCompressedTable.emplace_back(buffer.begin(), buffer.begin()+outputLen);
255  fCompressedTable.back().shrink_to_fit();
256 
257  compressedSize += outputLen;
258  }
259  std::cout << std::endl;
260  std::cout << "Shrunk " << N*sizeof(uint16_t)/(1024*1024) << "MB table to "
261  << compressedSize/(1024*1024) << "MB." << std::endl;
262  }
263 
264  //......................................................................
266  {
267  }
268 
269  //......................................................................
270  double ADCShapeFitTable::TNS(double tns0, int16_t adc1, int16_t adc2, int16_t adc3,
271  double& adcpeak, bool& goodTime, double& base) const
272  {
273  uint16_t offset;
274  double riseTime;
275  double fallTime ;
276  double preAmp;
277  const double sampleTime = fIsND ? 125 : 500; // ns
278 
279  if (fMode == 0){
280  riseTime = (fIsND ? 140 : 380)/sampleTime;
281  fallTime = (fIsND ? 4500 : 7000)/sampleTime;
282  preAmp = 0;
283 
284  if(!fIsND && !fIsMC) riseTime = 460/sampleTime;
285 
286  } else{
287  riseTime = (fIsND ? 116 : 424)/sampleTime;
288  fallTime = (fIsND ? 6729 : 10449)/sampleTime;
289  preAmp = ( fIsND ? 147000 : 110000)/sampleTime;
290 
291  //fd data 410 900
292  if(!fIsND && !fIsMC)
293  { riseTime = 410/sampleTime;
294  fallTime = 9000/sampleTime;
295  }
296  }
297  // These are the conditions for the best offset to be in the table
298  if(adc3 >= 0 && adc3 < kADC3Max &&
299  adc1 >= kADC1Min && adc1 < adc3 &&
300  adc2 >= adc3+kADC2RelMin && adc2 < adc3+kADC2RelMax){
301  ++fNHit;
302 
303  // And this is the offset it'll be at, based on aforementioned geometry
304  const int tableIdx = (kADC2RelMax-kADC2RelMin)* // Wedge thickness
305  (-kADC1Min*(adc3+1)+ // Rectangular part
306  (adc3*(adc3-1))/2+adc1)+ // Triangular part
307  (adc2-adc3-kADC2RelMin); // Distance through wedge
308 
309  assert(tableIdx >= 0 && tableIdx < int(kNumTableEntries));
310 
311  offset = Table(tableIdx);
312  goodTime = (fIsBad.count(tableIdx) == 0);
313 
314  // One time in 400, do the full calculation just to sanity check that the
315  // table (and our indexing into it) agrees. We have about a 5% miss rate
316  // anyway, so this isn't too much additional computation.
317 
318  if(adc3%20 == 0 && adc1%20 == 0){
319 
320  assert(ADCShapeFit(adc1, adc2, adc3, riseTime, fallTime, preAmp, goodTime, fMode) == offset);
321 
322  }
323  }
324  else{
325  // It's not in our table, just do the fit ourselves.
326  ++fNMiss;
327  offset = ADCShapeFit(adc1, adc2, adc3, riseTime, fallTime, preAmp, goodTime, fMode);
328  }
329 
330  // Empirically, if the best-fit t0 is before the first sample, then this is
331  // a bad fit. See docdb 12601.
332  if(offset < -kZeroOffsetSamples/kSamplesPerOffset) goodTime = false;
333 
334  // We need to extract the best normalization, so repeat the relevant part
335  // of the inner loop at that offset.
336  if(goodTime){
337  // The observations, in chronological order
338  const int16_t obs[kNumFineTimingADCPoints] = {0, adc1, adc2, adc3};
339  const double t0 = (double(offset)*kSamplesPerOffset+kZeroOffsetSamples);
340  double junk[kNumFineTimingADCPoints]; // Expectations
341  const double norm = GetExpectations(t0, riseTime, fallTime, preAmp, obs, junk, base, fMode);
342 
343  // The value the peak of the curve would have before scaling by norm
344  // Analytic calculation from the functional form in ReadoutSim
345  static const double peak = pow(riseTime/(riseTime+fallTime), riseTime/fallTime)-
346  pow(riseTime/(riseTime+fallTime), 1+riseTime/fallTime);
347 
348  adcpeak = norm*peak;
349  }
350  else{
351  // Unless the fit failed, in which case we're probably safest just doing
352  // this.
353  adcpeak = adc3;
354  return tns0;
355  }
356 
357  return tns0+(offset*kSamplesPerOffset+kZeroOffsetSamples)*sampleTime;
358  }
359 
360  //......................................................................
362  {
363  return fNHit ? fNHit/double(fNHit+fNMiss) : 0;
364  }
365 
366  //......................................................................
367  uint16_t ADCShapeFitTable::Table(int i) const
368  {
369  // If anything goes wrong and this compression/decompression doesn't
370  // reproduce the raw table, the checks in Calibrator::GetTNS(), where we
371  // occasionally re-fit a trace even though it's in the table, will catch
372  // it.
373 
374  static uint16_t buffer[kChunk];
375 
376  // Figure out which chunk 'i' is in
377  const int chunkIdx = i/kChunk;
378 
379  // Extract the relevant compressed chunk into our buffer
380  const unsigned int lenUsed =
381  LZ4_decompress_fast(&fCompressedTable[chunkIdx].front(),
382  (char*)buffer,
383  kChunk*sizeof(uint16_t));
384 
385  // Make sure it consumed the whole compressed chunk
386  assert(lenUsed == fCompressedTable[chunkIdx].size());
387 
388  // Return the element of the chunk corresponding to 'i'
389  return buffer[i%kChunk];
390  }
391 #endif // ADCSHAPEFIT_FUNC_ONLY
392 
393 } // namespace
std::vector< std::vector< char > > fCompressedTable
Chunks of LZ4-compressed data containing the lookup table.
Definition: ADCShapeFit.h:79
#define LZ4HC_CLEVEL_DEFAULT
Definition: lz4hc.h:48
constexpr T pow(T x)
Definition: pow.h:75
uint16_t Table(int i) const
Access the table entry at index i.
int LZ4_decompress_fast(const char *source, char *dest, int originalSize)
Definition: lz4.cxx:1261
static const unsigned int kNumTableEntries
Definition: ADCShapeFit.h:59
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
::xsd::cxx::tree::buffer< char > buffer
Definition: Database.h:179
static const int kADC2RelMax
Definition: ADCShapeFit.h:55
char int compressedSize
Definition: lz4.h:458
static const int kADC2RelMin
Definition: ADCShapeFit.h:54
Float_t Z
Definition: plot.C:38
CDPStorage service.
const double kZeroOffsetSamples
Definition: ADCShapeFit.cxx:50
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
double TNS(double tns0, int16_t adc1, int16_t adc2, int16_t adc3, double &adcpeak, bool &goodTime, double &base) const
static const int kADC1Min
Definition: ADCShapeFit.h:53
const unsigned int kNumFineTimingADCPoints
Definition: ADCShapeFit.cxx:45
uint16_t ADCShapeFit(int16_t adc1, int16_t adc2, int16_t adc3, double riseTime, double fallTime, double preAmp, bool &goodTime, int fMode)
OStream cout
Definition: OStream.cxx:6
ADCShapeFitTable(const std::string &fname, bool isND, bool isMC, int mode)
const int kChunk
Definition: ADCShapeFit.cxx:39
double TableHitFraction() const
What fraction of calls to TNS went via the table?
const double kSamplesPerOffset
Definition: ADCShapeFit.cxx:54
T sqr(T x)
Function to perform, and table to cache, timing fits to ADC values.
Definition: ADCShapeFit.cxx:23
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::unordered_set< int > fIsBad
Definition: ADCShapeFit.h:75
Float_t norm
int LZ4_compress_HC(const char *src, char *dst, int srcSize, int maxDstSize, int compressionLevel)
Definition: lz4hc.cxx:534
assert(nhit_max >=nhit_nbins)
int LZ4_compressBound(int isize)
Definition: lz4.cxx:395
double T
Definition: Xdiff_gwt.C:5
Float_t X
Definition: plot.C:38
#define W(x)
double GetExpectations(double t0, double riseTime, double fallTime, double preAmp, const int16_t *obs, double *exps, double &base, int mode)
Helper function for ADCShapeFit inner loop.
Definition: ADCShapeFit.cxx:64
const int kOffsetStep
Definition: ADCShapeFit.cxx:53
static const int kADC3Max
Definition: ADCShapeFit.h:52