AdcDPhase3x1x1LocalRoiBuilder_tool.cc
Go to the documentation of this file.
1 // AdcDPhase3x1x1LocalRoiBuilder_tool.cc
2 // christoph.alt@cern.ch
3 
4 //********************************
5 // This code is intended to be used on raw data. It determines ROI that should be excluded from pedestal and noise pattern calculation.
6 //
7 // Thode calculates a local pedestal in a sliding window and looks for ROI in the bins that follow this window.
8 // A ROI is defined when three criteria are met:
9 // 1. a minimum number of consecutive bins above a relatively low threhsold -> temporary ROI
10 // 2. inside temporary ROI: a minimum number of consecutive bins above a medium threshold
11 // 3. inside temporary ROI: at least one bin above a high threshold
12 // All thresholds are relative to the local pedestal of the sliding window. Threshold 1 and 2 can be expressed
13 // both in absolute ADC counts or as multiples of the standard deviation.
14 // The standard deviation is only calculated for the first window while the pedestal is re-calculated for each window.
15 //
16 //********************************
17 
19 #include <iostream>
20 #include <algorithm>
21 #include <sstream>
22 #include <iomanip>
24 
25 using std::vector;
26 using std::string;
27 using std::ostream;
28 using std::cout;
29 using std::endl;
30 using std::setw;
31 
32 //**********************************************************************
33 // Class methods.
34 //**********************************************************************
35 
37  m_LogLevel(ps.get<int>("LogLevel")),
38  m_BinsToAverageForPedestal(ps.get<AdcIndex>("BinsToAverageForPedestal")),
39  m_BinsToSkip(ps.get<AdcIndex>("BinsToSkip")),
40  m_UseStandardDeviation(ps.get<bool>("UseStandardDeviation")),
41  m_NConsecBinsAboveThreshold1(ps.get<AdcIndex>("NConsecBinsAboveThreshold1")),
42  m_NSigmaStart1(ps.get<AdcSignal>("NSigmaStart1")),
43  m_NSigmaEnd1(ps.get<AdcSignal>("NSigmaEnd1")),
44  m_NConsecBinsAboveThreshold2(ps.get<AdcIndex>("NConsecBinsAboveThreshold2")),
45  m_NSigmaStart2(ps.get<AdcSignal>("NSigmaStart2")),
46  m_NSigmaMax(ps.get<AdcSignal>("NSigmaMax")),
47  m_PadLow(ps.get<AdcIndex>("PadLow")),
48  m_PadHigh(ps.get<AdcIndex>("PadHigh")) {
49  const string myname = "AdcDPhase3x1x1LocalRoiBuilder::ctor: ";
50  if ( m_LogLevel > 0 ) {
51  cout << myname << " Configuration: " << endl;
52  cout << myname << " LogLevel: " << m_LogLevel << endl;
53  cout << myname << " BinsToAverageForPedestal: " << m_BinsToAverageForPedestal << endl;
54  cout << myname << " BinsToSkip: " << m_BinsToSkip << endl;
55  cout << myname << " UseStandardDeviation: " << m_UseStandardDeviation << endl;
56  cout << myname << " NConsecBinsAboveThreshold1: " << m_NConsecBinsAboveThreshold1 << endl;
57  cout << myname << " NSigmaStart1: " << m_NSigmaStart1 << endl;
58  cout << myname << " NSigmaEnd1: " << m_NSigmaEnd1 << endl;
59  cout << myname << " NConsecBinsAboveThreshold2: " << m_NConsecBinsAboveThreshold2 << endl;
60  cout << myname << " NSigmaStart2: " << m_NSigmaStart2 << endl;
61  cout << myname << " NSigmaMax: " << m_NSigmaMax << endl;
62  cout << myname << " PadLow: " << m_PadLow << endl;
63  cout << myname << " PadHigh: " << m_PadHigh << endl;
64  }
65 }
66 
67 //**********************************************************************
68 
69 
71 }
72 
73 //**********************************************************************
74 
76  const string myname = "AdcDPhase3x1x1LocalRoiBuilder:build: ";
77  if ( m_LogLevel >= 2 ) cout << myname << "Building ROIs for channel "
78  << data.channel() << "." << endl;
79  data.rois.clear();
80 
81  //Create dummy DataMap to return
82  DataMap res(0);
83  res.setInt("Test", 0);
84 
85  //Raw ADC.
86  AdcSignalVector sigs;
87  sigs = data.samples;
88 
89  // Build ROIS before padding and merging.
90  AdcFilterVector& signal = data.signal;
91  AdcRoiVector& rois = data.rois;
92  signal.clear();
93  signal.resize(sigs.size(), false);
94  bool inroi = false;
95  AdcSignal siglow1 = m_NSigmaEnd1;
96  AdcSignal sighigh1 = m_NSigmaStart1;
97  AdcSignal sighigh2 = m_NSigmaStart2;
98  AdcSignal sigmax = m_NSigmaMax;
99  AdcIndex nsig = sigs.size();
100 
101  int ROICount=0;
102  std::vector<AdcIndex> ROIStart;
103  std::vector<AdcIndex> ROIEnd;
104  ROIStart.clear();
105  ROIEnd.clear();
106 
107  if ( nsig < 1 ) {
108  if ( m_LogLevel >= 2 ) cout << myname << "Channel " << data.channel()
109  << " has no samples." << endl;
110  return res;
111  }
112 
113  //*****************************************
114  // FIRST ITERATION: BIN 0 TO LAST BIN ***
115  //*****************************************
116  double SumPedestal=0.;
117  double SumADCMinusPedestalSquared=0.;
118  double StandardDeviationPedestal=0.;
119  AdcIndex FirstEntryInPedestalSum=0;
120 
121  //Calculate pedestal for first m_BinsToAverageForPedestal ticks
123  {
124  SumPedestal+=sigs[isig];
125  }
126  FirstEntryInPedestalSum=m_BinsToSkip;
127 
128  //Calculate standard deviation for first m_BinsToAverageForPedestal ticks
130  {
131  for ( AdcIndex isig=m_BinsToSkip; isig<m_BinsToAverageForPedestal+m_BinsToSkip; ++isig )
132  {
133  SumADCMinusPedestalSquared+=pow(sigs[isig]-SumPedestal/m_BinsToAverageForPedestal,2);
134  }
135  StandardDeviationPedestal = sqrt(SumADCMinusPedestalSquared/(m_BinsToAverageForPedestal-1));
136  siglow1 = m_NSigmaEnd1*StandardDeviationPedestal;
137  sighigh1 = m_NSigmaStart1*StandardDeviationPedestal;
138  sighigh2 = m_NSigmaStart2*StandardDeviationPedestal;
139  }
140 
141  for ( AdcIndex isig = m_BinsToAverageForPedestal+m_BinsToSkip; isig<sigs.size(); ++isig )
142  {
143  AdcSignal sig = sigs[isig];
144  if ( inroi )
145  {
146  if ( sig > siglow1 + SumPedestal/m_BinsToAverageForPedestal && isig < sigs.size()-1)
147  {
148  signal[isig] = true;
149  }
150  else
151  {
152  ROIEnd.push_back(isig-1);
153  ROICount++;
154 
155  //check second criterion in this temporary ROI.
156  bool KeepThisROI = false;
157  AdcIndex m_NConsecBinsAboveThreshold2Temp = std::min(m_NConsecBinsAboveThreshold2,(AdcIndex)(ROIEnd[ROICount-1]-ROIStart[ROICount-1]+1));
158  AdcIndex NConsecBinsAboveThreshold2Count=0;
159 
160  for(AdcIndex isigroi = ROIStart[ROICount-1]; isigroi <= ROIEnd[ROICount-1]; isigroi++)
161  {
162  if( sigs[isigroi] >= sighigh2 + SumPedestal/m_BinsToAverageForPedestal )
163  {
164  NConsecBinsAboveThreshold2Count++;
165  if(NConsecBinsAboveThreshold2Count == m_NConsecBinsAboveThreshold2Temp)
166  {
167  KeepThisROI = true;
168  break;
169  }
170  }
171  else
172  {
173  NConsecBinsAboveThreshold2Count = 0;
174  }
175  }
176 
177  //if second criterion for this temporary ROI is met, check third criterion
178  if(KeepThisROI)
179  {
180  KeepThisROI = false;
181  for(AdcIndex isigroi = ROIStart[ROICount-1]; isigroi <= ROIEnd[ROICount-1]; isigroi++)
182  {
183  if( sigs[isigroi] >= sigmax + SumPedestal/m_BinsToAverageForPedestal)
184  {
185  KeepThisROI = true;
186  break;
187  }
188  }
189  }
190 
191  //if second or third criteria is not met, delete temporary ROI. Otherwise keep it.
192  if(!KeepThisROI)
193  {
194  for(AdcIndex isigroi = ROIStart[ROICount-1]; isigroi <= ROIEnd[ROICount-1]; isigroi++)
195  {
196  signal[isigroi] = false;
197  //recalculate pedesta, including the ones in the deleted ROI
198  SumPedestal -= sigs[FirstEntryInPedestalSum]; //remove last entry in pedestal sum
199  FirstEntryInPedestalSum++; //increase index for last entry in pedestal sum by 1
200  while(signal[FirstEntryInPedestalSum]) //check if the increased index was a signal. if yes, increase until index with no signal was found.
201  {
202  FirstEntryInPedestalSum++;
203  }
204  SumPedestal += sigs[isigroi]; //add current ADC count to pedestal sum
205  }
206 
207  ROIStart.pop_back();
208  ROIEnd.pop_back();
209  ROICount--;
210  }
211 
212  inroi = false;
213 
214  SumPedestal -= sigs[FirstEntryInPedestalSum];//remove last entry in pedestal sum
215 
216  FirstEntryInPedestalSum++; //increase index for last entry in pedestal sum by 1
217  while(signal[FirstEntryInPedestalSum]) //check if the increased index was a signal. if yes, increase until index with no signal was found.
218  {
219  FirstEntryInPedestalSum++;
220  }
221 
222  SumPedestal += sigs[isig]; //add current ADC count to sum
223  }
224  }
225  else
226  {
227  bool ROIStartIsHere = true;
228  for(AdcIndex isignext = isig; isignext < std::min((AdcIndex)sigs.size(),isig+m_NConsecBinsAboveThreshold1); isignext++)
229  {
230  if(sigs[isignext] < sighigh1 + SumPedestal/m_BinsToAverageForPedestal)
231  {
232  ROIStartIsHere = false;
233  break;
234  }
235  }
236 
237  if(ROIStartIsHere)
238  {
239  ROIStart.push_back(isig);
240  signal[isig] = true;
241  inroi = true;
242  }
243  else
244  {
245  SumPedestal -= sigs[FirstEntryInPedestalSum]; //remove last entry in pedestal sum
246 
247  FirstEntryInPedestalSum++; //increase index for last entry in pedestal sum by 1
248  while(signal[FirstEntryInPedestalSum]) //check if the increased index was a signal. if yes, increase until index with no signal was found.
249  {
250  FirstEntryInPedestalSum++;
251  }
252 
253  SumPedestal += sigs[isig]; //add current ADC count to sum
254  }
255  }
256  }
257 
258 //removal isolated signal at the end of the waveform
259 if(signal[sigs.size()-1] && !signal[sigs.size()-2])
260 {
261 signal[sigs.size()-1] = false;
262 }
263 
264  //******************************************
265  // SECOND ITERATION: LAST BIN TO BIN 0 ***
266  //******************************************
267  AdcIndex PedestalIndex=sigs.size()-1;
268  AdcIndex PedestalCounter=0;
269  inroi = false;
270  siglow1 = m_NSigmaEnd1;
271  sighigh1 = m_NSigmaStart1;
272  sighigh2 = m_NSigmaStart2;
273 
274  SumPedestal=0.;
275  SumADCMinusPedestalSquared=0.;
276  StandardDeviationPedestal=0.;
277 
278  FirstEntryInPedestalSum=0;
279 
280  ROIStart.clear();
281  ROIEnd.clear();
282  ROICount=0;
283 
284  //Calculate pedestal for last m_BinsToAverageForPedestal ticks
285  while(PedestalCounter < m_BinsToAverageForPedestal)
286  {
287  if(!signal[PedestalIndex])
288  {
289  if(PedestalCounter == 0) FirstEntryInPedestalSum=PedestalIndex; //remember first enntry of pedestal sum
290  SumPedestal+=sigs[PedestalIndex];
291  PedestalCounter++;
292  }
293  PedestalIndex--;
294  }
295 
296  //Calculate standard deviation for last m_BinsToAverageForPedestal ticks
298  {
299  PedestalCounter=0;
300  PedestalIndex=sigs.size()-1;
301  while(PedestalCounter < m_BinsToAverageForPedestal)
302  {
303  if(!signal[PedestalIndex])
304  {
305  SumADCMinusPedestalSquared+=pow(sigs[PedestalIndex]-SumPedestal/m_BinsToAverageForPedestal,2);
306  PedestalCounter++;
307  }
308  PedestalIndex--;
309  }
310  StandardDeviationPedestal = sqrt(SumADCMinusPedestalSquared/(m_BinsToAverageForPedestal-1));
311  siglow1 = m_NSigmaEnd1*StandardDeviationPedestal;
312  sighigh1 = m_NSigmaStart1*StandardDeviationPedestal;
313  sighigh2 = m_NSigmaStart2*StandardDeviationPedestal;
314  }
315 
316 
317  for ( AdcIndex isig = PedestalIndex; isig >= m_BinsToSkip && isig <= PedestalIndex; --isig )
318  {
319  if(signal[isig]) continue;
320  AdcSignal sig = sigs[isig];
321  if ( inroi )
322  {
323  if ( sig > siglow1 + SumPedestal/m_BinsToAverageForPedestal && isig > m_BinsToSkip)
324  {
325  signal[isig] = true;
326  }
327  else
328  {
329  ROIEnd.push_back(isig+1);
330  ROICount++;
331 
332  //check second criterion in this temporary ROI.
333  bool KeepThisROI = false;
334  AdcIndex m_NConsecBinsAboveThreshold2Temp = std::min(m_NConsecBinsAboveThreshold2,(AdcIndex)(ROIEnd[ROICount-1]-ROIStart[ROICount-1]+1));
335  AdcIndex NConsecBinsAboveThreshold2Count=0;
336 
337  for(AdcIndex isigroi = ROIStart[ROICount-1]; isigroi >= ROIEnd[ROICount-1]; isigroi--)
338  {
339  if( sigs[isigroi] >= sighigh2 + SumPedestal/m_BinsToAverageForPedestal )
340  {
341  NConsecBinsAboveThreshold2Count++;
342  if(NConsecBinsAboveThreshold2Count == m_NConsecBinsAboveThreshold2Temp)
343  {
344  KeepThisROI = true;
345  break;
346  }
347  }
348  else
349  {
350  NConsecBinsAboveThreshold2Count = 0;
351  }
352  }
353 
354  //if second criterion for this temporary ROI is met, check third criterion
355  if(KeepThisROI)
356  {
357  KeepThisROI = false;
358  for(AdcIndex isigroi = ROIStart[ROICount-1]; isigroi >= ROIEnd[ROICount-1]; isigroi--)
359  {
360  if( sigs[isigroi] >= sigmax + SumPedestal/m_BinsToAverageForPedestal )
361  {
362  KeepThisROI = true;
363  break;
364  }
365  }
366  }
367 
368  //if second or third criteria is not met, delete temporary ROI. Otherwise keep it.
369  if(!KeepThisROI)
370  {
371  for(AdcIndex isigroi = ROIStart[ROICount-1]; isigroi >= ROIEnd[ROICount-1]; isigroi--)
372  {
373  signal[isigroi] = false;
374  //recalculate pedestal, including the ones in the deleted ROI
375  SumPedestal -= sigs[FirstEntryInPedestalSum]; //remove last entry in pedestal sum
376  FirstEntryInPedestalSum--; //decreas index for last entry in pedestal sum by 1
377  while(signal[FirstEntryInPedestalSum]) //check if the decreased index was a signal. if yes, increase until index with no signal was found.
378  {
379  FirstEntryInPedestalSum--;
380  }
381  SumPedestal += sigs[isigroi]; //add current ADC count to sum
382  }
383 
384  ROIStart.pop_back();
385  ROIEnd.pop_back();
386  ROICount--;
387  }
388 
389  inroi = false;
390  SumPedestal -= sigs[FirstEntryInPedestalSum];//remove last entry in pedestal sum
391 
392  FirstEntryInPedestalSum--; //increase index for last entry in pedestal sum by 1
393  while(signal[FirstEntryInPedestalSum]) //check if the increased index was a signal. if yes, increase until index with no signal was found.
394  {
395  FirstEntryInPedestalSum--;
396  }
397 
398  SumPedestal += sigs[isig]; //add current ADC count to sum
399 
400  }
401  }
402  else
403  {
404  bool ROIStartIsHere = true;
405  for(AdcIndex isignext = isig; isignext >= (AdcIndex)std::max((int)m_BinsToSkip,(int)isig-(int)m_NConsecBinsAboveThreshold1+1) && isignext <= PedestalIndex ; isignext--)
406  {
407  if(sigs[isignext] < sighigh1 + SumPedestal/m_BinsToAverageForPedestal)
408  {
409  ROIStartIsHere = false;
410  break;
411  }
412  }
413 
414  if(ROIStartIsHere)
415  {
416  ROIStart.push_back(isig);
417  signal[isig] = true;
418  inroi = true;
419  }
420  else
421  {
422  SumPedestal -= sigs[FirstEntryInPedestalSum]; //remove first entry in pedestal sum
423  FirstEntryInPedestalSum--; //decrease index for last entry in pedestal sum by 1
424 
425  while(signal[FirstEntryInPedestalSum]) //check if the increased index was a signal. if yes, increase until index with no signal was found.
426  {
427  FirstEntryInPedestalSum--;
428  }
429  SumPedestal += sigs[isig]; //add current ADC count to sum
430  }
431  }
432  }
433 
434 //removal isolated signal at the beginning of the waveform
435 if(signal[m_BinsToSkip] && !signal[m_BinsToSkip+1])
436 {
437 signal[m_BinsToSkip] = false;
438 }
439 
440  // Fill the unpadded ROIs.
441  data.roisFromSignal();
442  // Display ROIs before padding and merging.
443  if ( m_LogLevel >= 3 ) {
444  cout << myname << " ROIs before merge (size = " << rois.size() << "):" << endl;
445  for ( const AdcRoi& roi : rois ) {
446  cout << myname << setw(8) << roi.first << " " << setw(8) << roi.second << endl;
447  }
448  } else if ( m_LogLevel >= 2 ) {
449  cout << myname << " ROI count before merge: " << data.rois.size() << endl;
450  }
451  if ( rois.size() == 0 ) return res;
452  // Pad ROIs.
453  unsigned int isig1 = 0;
454  unsigned int isig2 = 0;
455  for ( AdcRoi roi : rois ) {
456  isig2 = roi.first;
457  isig1 = isig2 > m_PadLow ? isig2 - m_PadLow : 0;
458  for ( unsigned int isig=isig1; isig<isig2; ++isig ) signal[isig] = true;
459  isig1 = roi.second + 1;
460  isig2 = isig1 + m_PadHigh;
461  if ( isig2 > nsig ) isig2 = nsig;
462  for ( unsigned int isig=isig1; isig<isig2; ++isig ) signal[isig] = true;
463  }
464  // Fill the final ROIs.
465  data.roisFromSignal();
466  // Display final ROIs.
467  if ( m_LogLevel >= 3 ) {
468  cout << myname << " ROIs after merge (size = " << rois.size() << "):" << endl;
469  for ( const AdcRoi& roi : rois ) {
470  cout << myname << setw(8) << roi.first << " " << setw(8) << roi.second << endl;
471  }
472  } else if ( m_LogLevel >= 2 ) {
473  cout << myname << " ROI count after merge: " << data.rois.size() << endl;
474  }
475  return res;
476 }
477 //**********************************************************************
478 
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
std::string string
Definition: nybbler.cc:12
float AdcSignal
Definition: AdcTypes.h:21
constexpr T pow(T x)
Definition: pow.h:72
struct vector vector
std::pair< AdcIndex, AdcIndex > AdcRoi
Definition: AdcTypes.h:54
AdcRoiVector rois
void setInt(Name name, int val)
Definition: DataMap.h:131
static int max(int a, int b)
static constexpr double ps
Definition: Units.h:99
unsigned int AdcIndex
Definition: AdcTypes.h:15
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
std::vector< AdcRoi > AdcRoiVector
Definition: AdcTypes.h:55
Channel channel() const
AdcFilterVector signal
std::vector< bool > AdcFilterVector
Definition: AdcTypes.h:27
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
std::vector< AdcSignal > AdcSignalVector
Definition: AdcTypes.h:22
AdcDPhase3x1x1LocalRoiBuilder(fhicl::ParameterSet const &ps)
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
int bool
Definition: qglobal.h:345
DataMap update(AdcChannelData &acd) const override
AdcSignalVector samples
QTextStream & endl(QTextStream &s)