FPGAAlgorithms.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////////////
2 /// \brief TODO
3 /// \author ??
4 /// \date
5 /////////////////////////////////////////////////////////////////////////////////////
6 
8 
11 #include "cetlib_except/exception.h"
12 #include "cetlib/search_path.h"
16 
17 #include "DAQDataFormats/NanoSliceVersionConvention.h"
18 #include "Geometry/Geometry.h"
19 #include "NovaDAQConventions/DAQConventions.h"
22 
23 #include "sys/stat.h"
24 
25 #include <algorithm>
26 #include <limits>
27 
28 #include "TSystem.h"
29 #include "TFile.h"
30 
31 namespace rsim
32 {
34  CommonParameters(pset)
35  ,fThresholdADC(0)
36  ,fThresholdPDF(0)
37  ,fThresholds(0)
38  ,fPlane(0)
39  ,fCell(0)
40  {
41  //cannot get thresholds from the database and randomly
43  throw cet::exception("Configuration") << "FPGAAlgorithm: UseRandomThreshold and UseDBThreshold cannot both be true \n";
44  }
45 
46  //fetch the threshold PDF if we are setup to use random thresholds
48  // constructor decides if initialized value is a path or an environment variable
49  cet::search_path sp("FW_SEARCH_PATH");
50 
51  sp.find_file(pset.get<std::string>("ThresholdFile"), fThresholdFile);
52 
53  struct stat sb;
54  if (fThresholdFile.empty() || stat(fThresholdFile.c_str(), &sb) != 0){
55  // failed to resolve the file name
56  throw cet::exception("FileNotOpen") << "FPGAAlgorithms: Unable to find " << fThresholdFile << " with the threshold PDF \n "
57  <<__FILE__ << ":" << __LINE__ << "\n";
58  }
59 
60  if (gSystem->AccessPathName(fThresholdFile.c_str())){
61  throw cet::exception("FileNotOpen") << "FPGAAlgorithms: Unable to open " << fThresholdFile << " to get the threshold PDF \n";
62  }
63 
64  TFile f(fThresholdFile.c_str());
65  TH1D *htemp = (TH1D*)f.Get("thresholds");
66  if (!htemp){
67  throw cet::exception("FatalRootError") << "FPGAAlgorithms: Unable to get thresholds histogram from " << fThresholdFile << " \n";
68  }
69  fThresholdPDF = (TH1D*)htemp->Clone("ThresholdPDF");
70  fThresholdPDF->SetDirectory(0);
71  f.Close();
72  }
73  }
74 
76  {
77  if (fThresholdPDF) delete fThresholdPDF;
78  }
79 
82  IFPGAAlgorithm(pset),
83  fLookbackSamples (pset.get<int >("FPGA_DCS_LookbackSamples")),
84  fRetriggerThreshold(pset.get<double>("FPGA_DCS_RetriggerThreshold")),
85  fVersion (pset.get<int >("FPGA_DCS_NanosliceVersion"))
86  {
87  fThresholdADC = pset.get<double>("FPGA_DCS_ThresholdADC");
88  }
89 
92  IFPGAAlgorithm(pset),
93  fBaselineTime(pset.get<double>("FPGA_MatchedFiltering_BaselineTime")),
94  fEndTime (pset.get<double>("FPGA_MatchedFiltering_EndTime"))
95  {
96  fThresholdADC = pset.get<double>("FPGA_MatchedFiltering_ThresholdMatchValue");
97  }
98 
100  IFPGAAlgorithm(pset)
101  {
102  fUseDBThresholds = false;
103  fThresholdADC = 0;
104  }
105 
106 
108  {
109  //should correctly not use the database for DumpASICOutput
111 
112  return fThresholdMap[fPlane][fCell];
113  }
114 
116  {
117  fThresholdMap.clear();
118 
120 
122  // get the RandomNumberGenerator service
124  CLHEP::HepRandomEngine &engine = rng->getEngine();
125  RandHisto randHisto(engine);
126  double thresh(0.0), dummy1(0.0), dummy2(0.0);
127 
128  for (unsigned int plane = 0; plane < geom->NPlanes(); plane++){
129  std::vector< int > plane_thresholds;
130  for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++){
131  randHisto.GetRandom( fThresholdPDF, thresh, dummy1, dummy2 );
132  plane_thresholds.push_back(thresh);
133  }
134  fThresholdMap.push_back(plane_thresholds);
135  }
136  }
137  else if (fUseDBThresholds){
139  for (unsigned int plane = 0; plane < geom->NPlanes(); plane++) {
140  std::vector< int > plane_thresholds;
141  for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++){
142  bool ok(false);
143  geo::OfflineChan oChan(plane, cell);
144  //int threshold = RH->GetThreshold(oChan, ok);
145  int threshold = RH->GetPedestal(oChan, ok);
146  //GetPedestal does not set the ok bool
147  if (threshold >= 10) plane_thresholds.push_back( threshold );
148  else plane_thresholds.push_back(4095);
149  }
150  fThresholdMap.push_back(plane_thresholds);
151  }
152  }
153 
154  if (!fThresholds){
156  fThresholds = tfs->make<TH1D>("fThresholds", ";Thresholds", 501, 0, 500);
157  for (unsigned int i = 0; i < fThresholdMap.size(); ++i){
158  for (unsigned int j = 0; j < fThresholdMap[i].size(); ++j){
159  fThresholds->Fill( fThresholdMap[i][j] );
160  }
161  }
162  }
163  return;
164  }
165 
166 
168  {
169  const std::string algo = pset.get<std::string>("FPGAAlgorithm");
170 
171 
172  if(algo == "DualCorrelatedSampling")
173  return new FPGA_DualCorrelatedSampling(pset);
174  if(algo == "FPGA_MatchedFiltering")
175  return new FPGA_MatchedFiltering(pset);
176  if(algo == "DumpASICOutput")
177  return new FPGA_DumpASICOutput(pset);
178 
179 
180  assert(0 && "Unknown FPGA algorithm");
181 
182  return nullptr;
183  }
184 
185 
186  //-------------- DualCorrelatedSampling --------------
187  std::vector<rawdata::RawDigit> FPGA_DualCorrelatedSampling::
188  ASICToDigits(double *ASIC_Output,
189  int firstDigitizationClocktick)
190  {
191  std::vector<rawdata::RawDigit> digitlist(0);
192 
193  const double threshold = ThresholdADC();
194 
196 
197  const int nSamples = conv.getNSamples(fVersion);
198  const int nPretrig = conv.getNPretriggeredSamples(fVersion);
199 
200  // Fill DCS space with the differences, skipping over one each time:
201  // Digitizations: Z0 Z1 Z2 ...
202  // DCS space: Z2-Z0 Z3-Z1 Z4-Z2 ...
203  //
204  // Look through DCS space searching for above-threshold values. When you
205  // find one, find the biggest DCS value that occurs and return that as a
206  // hit (with the time from the first >thresh point and the amplitude from
207  // the peak.)
208 
209  bool triggered = false;
210  double peakDCSValue = 0;
211  int peakDCSTime = 0;
212 
213  double DCS_Space[fNumClockticksInSpill];
214 
215  // Make sure neither the lookback (which actually looks forwards from the
216  // baseline to the peak the way we have it set up) nor the pretrig or any
217  // postrig samples can overflow the array bounds. This does cause us to
218  // miss hits right at the edge of the trigger window, which doesn't happen
219  // in real data. Should be a sub-percent effect on the overall rate, at the
220  // edges where hits can be legitimately cut off too.
221  const int iMin = firstDigitizationClocktick + std::max(nPretrig-fLookbackSamples, 0)*fNumClockticksPerDigitization;
222  const int iMax = fNumClockticksInSpill - (nSamples-nPretrig+fLookbackSamples-1)*fNumClockticksPerDigitization;
223 
224  for (int iClocktick = iMin; iClocktick < iMax; iClocktick += fNumClockticksPerDigitization) {
225 
226  DCS_Space[iClocktick] = ASIC_Output[iClocktick + fLookbackSamples*fNumClockticksPerDigitization] - ASIC_Output[iClocktick];
227 
228  // are we triggered?
229  double value = DCS_Space[iClocktick];
230 
231  if (value > threshold) {
232  triggered = true;
233  }
234  else if (triggered && value < fRetriggerThreshold*threshold) {
235  // End of triggered period declared when back below retrigger
236  // threshold. Report what we found and reset.
237 
238  rawdata::RawDigit digit;
239  digit.SetVersion(fVersion);
240 
241  digit.SetTDC((peakDCSTime+fLookbackSamples*fNumClockticksPerDigitization)<<2); // shift by two bits to match real DAQ
242 
243  if(fVersion == 0){
244  // Version zero is a single DCS value
245  digit.SetADC(0, peakDCSValue);
246  }
247  else{
248  // Other version are multiple raw values
249  const int dt = fNumClockticksPerDigitization;
250  for(int i = fLookbackSamples-nPretrig; i < nSamples+fLookbackSamples-nPretrig; ++i){
251  const int idx = peakDCSTime+i*dt;
252  // Check we got iMin and iMax correct above
253  assert(idx >= 0 && idx < fNumClockticksInSpill);
254  digit.SetADC(digit.NADC(), ASIC_Output[idx]);
255  }
256  }
257 
258  digitlist.push_back(digit);
259 
260  triggered = false;
261  peakDCSValue = 0;
262  peakDCSTime = 0;
263  }
264 
265  // looking for peak?
266  if (triggered) {
267  if (value > peakDCSValue) {
268  peakDCSTime = iClocktick;
269  peakDCSValue = value;
270  }
271  }
272 
273  }
274 
275  return digitlist;
276  }
277 
278 
279  //-------------- MatchedFiltering --------------
280  std::vector<rawdata::RawDigit> FPGA_MatchedFiltering::
281  ASICToDigits(double *ASIC_Output,
282  int firstDigitizationClocktick)
283  {
284  std::vector<rawdata::RawDigit> digitlist(0);
285  rawdata::RawDigit digit;
286 
287  const double threshold = ThresholdADC();
288 
289  // Convert the algo parameters into clocktick units (rounded to the nearest
290  // numClockticksPerDigitization) or into ADC counts per DCS gap (from PE/ns)
291  //
292  // The clocktick numbers here;
293  // BaselineWindow: Number of clockticks for baseline averaging
294  // MatchingWindow: Number of clockticks for matching calculation
295 
296  double MatchStartTime = 0; // Since f(t)=0 for t<0, no reason to include t<0 in sum
297  int MatchStartTimeFineticks = int(MatchStartTime / fClocktick);
298 
299  int BaselineWindow =
301  std::max( 1, int( (MatchStartTime - fBaselineTime)
303  int MatchingWindow =
305  std::max( 2, int( (fEndTime - MatchStartTime)
307 
308 
309  // Matched filtering provides precise time determination, but it doesn't
310  // help with the magnitude of the hit. We'll need to use the ASIC output
311  // directly to get that.
312  //
313  // Loop over the relevant digitizations for this channel
314 
315  bool triggered = false;
316  double peakMatchValue = 0;
317  int peakMatchTime = 0;
318 
319  // RBP debug bool dumpit = 0;
320  // RBP debug double tmpvalue[fNumClockticksInSpill];
321  // RBP debug for (int ii=0;ii<fNumClockticksInSpill;ii++) tmpvalue[ii] = 0;
322 
323  for (int iClocktick = firstDigitizationClocktick;
324  iClocktick < fNumClockticksInSpill - BaselineWindow - MatchingWindow; // ensure there's room left in the window (maybe change this later)
325  iClocktick += fNumClockticksPerDigitization) {
326 
327  // The first BaselineWindow ticks are devoted to calculating the baseline
328  // average
329  double baseline = 0;
330  int baseline_count = 0;
331  for (int iBaselinetick = 0;
332  iBaselinetick < BaselineWindow;
333  iBaselinetick += fNumClockticksPerDigitization) {
334  baseline += ASIC_Output[iClocktick+iBaselinetick];
335  baseline_count++;
336  }
337  baseline /= baseline_count; // With what precision can the FPGA do this?
338  // Also, not checking for div-by-0 (shouldn't be possible).
339 
340  // Loop through possible relative alignments of the test signal.
341  // There is no real reason this has to have anything to do with
342  // the clock. This finer looping can be on any scale. But, this
343  // is how it was described to me.
344  for (int iFinetick = 0;
345  iFinetick < fNumClockticksPerDigitization;
346  iFinetick++) {
347 
348  // Calculate the match value (inner product). Use the
349  // PerfectASICCurve function, although the ASIC output
350  // generation and this FPGA algorithm would, in general,
351  // use different functions (since the *real* ASIC could
352  // behave non-ideally.)
353  double value = 0;
354  double norm1 = 0;
355  double norm2 = 0;
356  int count = 0;
357  for (int iM=0; iM<MatchingWindow; iM+=fNumClockticksPerDigitization) {
358  double y1 = (ASIC_Output[iClocktick + BaselineWindow + iM] - baseline);
359  double y2 = PerfectASICCurve(fClocktick * (iM - iFinetick) + MatchStartTime);
360  value += y1 * y2;
361  norm1 += y1 * y1;
362  norm2 += y2 * y2;
363  count++;
364  }
365 
366  // Remove the effect of different numbers of digitizations used
367  // and simultaneously make the value represent a simple scaling of
368  // the input:
369  value /= norm2;
370 
371  // RBP debug int tmptime = iClocktick+BaselineWindow+iFinetick-MatchStartTimeFineticks;
372  // RBP debug if (tmptime<fNumClockticksInSpill&&tmptime>=0) {
373  // RBP debug tmpvalue[tmptime] = value;
374  // RBP debug }
375  // RBP debug else {
376  // RBP debug cout << "RBP - tmptime = " << tmptime << endl;
377  // RBP debug }
378 
379  // are we triggered?
380  if (value > threshold) {
381  triggered = true;
382  }
383  else if (triggered && value < 0.67*threshold) {
384  // End of triggered period declared when back below two-thirds
385  // threshold. Report what we found and reset. The "0" says to put
386  // the ADC/TDC info into the first element of the RawDigit vector.
387  // The vector nature of ADC and TDC in RawDigit isn't used at the
388  // moment. It could be used to store waveforms.
389  // RBP debug if (1) {
390  // RBP debug std::cout << "==== " << peakMatchTime << endl; RBP debug
391  // dumpit = 1; RBP debug }
392  digit.SetADC(0, uint16_t(peakMatchValue));
393  digit.SetTDC((unsigned int)peakMatchTime<<2); // shift by two bits to match real DAQ
394  digitlist.push_back(digit);
395  triggered = false;
396  peakMatchValue = 0;
397  peakMatchTime = 0;
398  }
399 
400  // looking for peak?
401  if (triggered) {
402  if (value > peakMatchValue) {
403  // For time, subtract MatchStartTimeFineticks so that
404  // MatchStartTime doesn't appear as a global shift. (Wouldn't
405  // really matter, but might be less confusing this way.) Also add
406  // in the baseline window to get the peak more or less in line with
407  // the ASIC output's ramp up.
408  peakMatchTime = iClocktick+BaselineWindow+iFinetick-MatchStartTimeFineticks;
409  peakMatchValue = value;
410  }
411  }
412 
413  }
414  }
415 
416  // RBP debug // RBP
417  // RBP debug if (dumpit&&0) {
418  // RBP debug FPGA_DumpASICOutput(ASIC_Output);
419  // RBP debug cout << "++" << endl;
420  // RBP debug FPGA_DumpASICOutput(&tmpvalue[0]);
421  // RBP debug }
422 
423 
424  return digitlist;
425  }
426 
428  {
429  // TODO: This is copy-and-pasted from SimpleReadout, which sucks
430 
431  double R = fASICRiseTime;
432  double F = fASICFallTime;
433 
434  static double peak = 0;
435  if(peak == 0) {
436  peak = pow(R/(R+F),R/F)-pow(R/(R+F),1+R/F);
437  }
438 
439  if (t<0) {
440  return 0;
441  }
442  else {
443  return exp(-t/F)*(1-exp(-t/R))/peak;
444  }
445  }
446 
447 
448  //-------------- DumpASICOutput (for debug) --------------
449  std::vector<rawdata::RawDigit> FPGA_DumpASICOutput::
450  ASICToDigits(double *ASIC_Output,
451  int firstDigitizationClocktick)
452  {
453  std::vector<rawdata::RawDigit> digitlist(0);
454  for (int iClocktick=0; iClocktick<fNumClockticksInSpill; iClocktick++) {
455  double tickTime = iClocktick*fClocktick;
456  mf::LogDebug("SimpleReadout") << tickTime << " " << ASIC_Output[iClocktick];
457  }
458 
459  return digitlist;
460  }
461 
462 } // end namespace
T max(const caf::Proxy< T > &a, T b)
#define F(x, y, z)
int fNumClockticksPerDigitization
How many ticks between each ADC.
double fASICRiseTime
ASIC shaping curve: rise time, ns.
friend IFPGAAlgorithm * CreateFPGAAlgorithm(const fhicl::ParameterSet &)
Create and configure the correct algorithm based on pset.
virtual std::vector< rawdata::RawDigit > ASICToDigits(double *ASIC_Output, int offset)
Float_t y1[n_points_granero]
Definition: compare.C:5
void SetTDC(int32_t iTDC)
Definition: RawDigit.cxx:116
constexpr T pow(T x)
Definition: pow.h:75
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
IFPGAAlgorithm(const fhicl::ParameterSet &pset)
Use CreateFPGAAlgorithm to create.
const PlaneGeo * Plane(unsigned int i) const
std::string find_file(std::string const &filename) const
void SetVersion(uint8_t v)
Definition: RawDigit.h:101
Common configuration params for SimpleReadout, FPGAAlgorithms, NoiseMaker.
Definition: Cand.cxx:23
void SetADC(uint32_t i, int16_t iADC)
Definition: RawDigit.cxx:108
base_engine_t & getEngine() const
FPGA_DualCorrelatedSampling(const fhicl::ParameterSet &pset)
Use CreateFPGAAlgorithm to create.
uint32_t getNSamples(const version_t ver) const
Get the number of samples.
int GetPedestal(const geo::OfflineChan &, bool &)
Definition: RunHistory.cxx:632
uint32_t getNPretriggeredSamples(const version_t ver) const
Get number of pretriggered samples.
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual std::vector< rawdata::RawDigit > ASICToDigits(double *ASIC_Output, int offset)
bool fUseDBThresholds
Use thresholds from DB instead of fcl defaults.
T get(std::string const &key) const
Definition: ParameterSet.h:231
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Interface implemented by all FPGA algorithms.
std::string fThresholdFile
File containing thresholds for the random threshold option.
#define R(x)
const double j
Definition: BetheBloch.cxx:29
double ThresholdADC() const
bool fUseRandomThresholds
Use thresholds drawn from a histogram instead of fcl defaults.
virtual std::vector< rawdata::RawDigit > ASICToDigits(double *ASIC_Output, int offset)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double PerfectASICCurve(double t) const
uint32_t NADC() const
Definition: RawDigit.h:81
std::vector< std::vector< int > > fThresholdMap
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
A (plane, cell) pair.
Definition: OfflineChan.h:17
Can be used as either a member holding configurations, or a mix-in.
void geom(int which=0)
Definition: geom.C:163
int fNumClockticksInSpill
number of clockticks in a spill
assert(nhit_max >=nhit_nbins)
void GetRandom(TH1 *histo, double &x, double &y, double &z)
Definition: RandHisto.cxx:8
FPGA_DumpASICOutput(const fhicl::ParameterSet &pset)
Use CreateFPGAAlgorithm to create.
unsigned int NPlanes() const
FPGA_MatchedFiltering(const fhicl::ParameterSet &pset)
Use CreateFPGAAlgorithm to create.
Encapsulate the geometry of one entire detector (near, far, ndos)
double fClocktick
digitization clock period, ns
double fASICFallTime
ASIC shaping curve: fall time, ns.