CMS 3D CMS Logo

SiStripFedZeroSuppression.cc
Go to the documentation of this file.
2 
7 
8 //#define DEBUG_SiStripZeroSuppression_
9 //#define ML_DEBUG
10 using namespace std;
11 
13  uint32_t n_cache_id = es.get<SiStripNoisesRcd>().cacheIdentifier();
14  uint32_t t_cache_id = es.get<SiStripThresholdRcd>().cacheIdentifier();
15 
16  if (n_cache_id != noise_cache_id) {
17  es.get<SiStripNoisesRcd>().get(noiseHandle);
18  noise_cache_id = n_cache_id;
19  }
20  if (t_cache_id != threshold_cache_id) {
21  es.get<SiStripThresholdRcd>().get(thresholdHandle);
22  threshold_cache_id = t_cache_id;
23  }
24 }
25 
26 void SiStripFedZeroSuppression::suppress(const std::vector<SiStripDigi>& in,
27  std::vector<SiStripDigi>& selectedSignal,
28  uint32_t detID) {
29  suppress(in, selectedSignal, detID, noiseHandle, thresholdHandle);
30 }
31 
32 void SiStripFedZeroSuppression::suppress(const std::vector<SiStripDigi>& in,
33  std::vector<SiStripDigi>& selectedSignal,
34  uint32_t detID,
35  edm::ESHandle<SiStripNoises>& noiseHandle,
36  edm::ESHandle<SiStripThreshold>& thresholdHandle) {
37  int i;
38  int inSize = in.size();
39  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
40  SiStripThreshold::Range detThRange = thresholdHandle->getRange(detID);
41 
42  // reserving more than needed, but quicker than one at a time
43  selectedSignal.clear();
44  selectedSignal.reserve(inSize);
45  for (i = 0; i < inSize; i++) {
46  //Find adc values for neighbouring strips
47  const uint32_t strip = (uint32_t)in[i].strip();
48 
49  adc = in[i].adc();
50 
51  SiStripThreshold::Data thresholds = thresholdHandle->getData(strip, detThRange);
52  theFEDlowThresh = static_cast<int16_t>(thresholds.getLth() * noiseHandle->getNoiseFast(strip, detNoiseRange) + 0.5);
53  theFEDhighThresh =
54  static_cast<int16_t>(thresholds.getHth() * noiseHandle->getNoiseFast(strip, detNoiseRange) + 0.5);
55 
56  adcPrev = -9999;
57  adcNext = -9999;
58  adcPrev2 = -9999;
59  adcNext2 = -9999;
60 
61  /*
62  Since we are not running on
63  Raw data we need to initialize
64  all the FED threshold
65  */
66 
67  theNextFEDlowThresh = theFEDlowThresh;
68  theNext2FEDlowThresh = theFEDlowThresh;
69  thePrevFEDlowThresh = theFEDlowThresh;
70  thePrev2FEDlowThresh = theFEDlowThresh;
71  theNeighFEDlowThresh = theFEDlowThresh;
72 
73  theNextFEDhighThresh = theFEDhighThresh;
74  thePrevFEDhighThresh = theFEDhighThresh;
75  theNeighFEDhighThresh = theFEDhighThresh;
76 
77  if (((strip) % 128) == 127) {
78  adcNext = 0;
79  theNextFEDlowThresh = 9999;
80  theNextFEDhighThresh = 9999;
81  } else if (i + 1 < inSize && in[i + 1].strip() == strip + 1) {
82  adcNext = in[i + 1].adc();
83  SiStripThreshold::Data thresholds_1 = thresholdHandle->getData(strip + 1, detThRange);
84  theNextFEDlowThresh =
85  static_cast<int16_t>(thresholds_1.getLth() * noiseHandle->getNoiseFast(strip + 1, detNoiseRange) + 0.5);
86  theNextFEDhighThresh =
87  static_cast<int16_t>(thresholds_1.getHth() * noiseHandle->getNoiseFast(strip + 1, detNoiseRange) + 0.5);
88  if (((strip) % 128) == 126) {
89  adcNext2 = 0;
90  theNext2FEDlowThresh = 9999;
91  } else if (i + 2 < inSize && in[i + 2].strip() == strip + 2) {
92  adcNext2 = in[i + 2].adc();
93  theNext2FEDlowThresh = static_cast<int16_t>(thresholdHandle->getData(strip + 2, detThRange).getLth() *
94  noiseHandle->getNoiseFast(strip + 2, detNoiseRange) +
95  0.5);
96  }
97  }
98 
99  if (((strip) % 128) == 0) {
100  adcPrev = 0;
101  thePrevFEDlowThresh = 9999;
102  thePrevFEDhighThresh = 9999;
103  } else if (i - 1 >= 0 && in[i - 1].strip() == strip - 1) {
104  adcPrev = in[i - 1].adc();
105  SiStripThreshold::Data thresholds_1 = thresholdHandle->getData(strip - 1, detThRange);
106  thePrevFEDlowThresh =
107  static_cast<int16_t>(thresholds_1.getLth() * noiseHandle->getNoiseFast(strip - 1, detNoiseRange) + 0.5);
108  thePrevFEDhighThresh =
109  static_cast<int16_t>(thresholds_1.getHth() * noiseHandle->getNoiseFast(strip - 1, detNoiseRange) + 0.5);
110  if (((strip) % 128) == 1) {
111  adcPrev2 = 0;
112  thePrev2FEDlowThresh = 9999;
113  } else if (i - 2 >= 0 && in[i - 2].strip() == strip - 2) {
114  adcPrev2 = in[i - 2].adc();
115  thePrev2FEDlowThresh = static_cast<int16_t>(thresholdHandle->getData(strip - 2, detThRange).getLth() *
116  noiseHandle->getNoiseFast(strip - 2, detNoiseRange) +
117  0.5);
118  }
119  }
120 
121  if (adcNext <= adcPrev) {
122  adcMaxNeigh = adcPrev;
123  theNeighFEDlowThresh = thePrevFEDlowThresh;
124  theNeighFEDhighThresh = thePrevFEDhighThresh;
125  } else {
126  adcMaxNeigh = adcNext;
127  theNeighFEDlowThresh = theNextFEDlowThresh;
128  theNeighFEDhighThresh = theNextFEDhighThresh;
129  }
130 
131  if (isAValidDigi()) {
132  selectedSignal.push_back(SiStripDigi(strip, adc));
133  }
134  }
135 }
136 
138  const uint32_t detID = out.id;
139  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
140  SiStripThreshold::Range detThRange = thresholdHandle->getRange(detID);
141 #ifdef DEBUG_SiStripZeroSuppression_
142  if (edm::isDebugEnabled())
143  LogTrace("SiStripZeroSuppression")
144  << "[SiStripFedZeroSuppression::suppress] Zero suppression on edm::DetSet<SiStripRawDigi>: detID " << detID
145  << " size = " << in.data.size();
146 #endif
147  edm::DetSet<SiStripRawDigi>::const_iterator in_iter = in.data.begin();
148  for (; in_iter != in.data.end(); in_iter++) {
149  const uint32_t strip = (uint32_t)(in_iter - in.data.begin());
150 
151 #ifdef DEBUG_SiStripZeroSuppression_
152  if (edm::isDebugEnabled())
153  LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] detID= " << detID
154  << " strip= " << strip << " adc= " << in_iter->adc();
155 #endif
156  adc = in_iter->adc();
157 
158  SiStripThreshold::Data thresholds = thresholdHandle->getData(strip, detThRange);
159  theFEDlowThresh = static_cast<int16_t>(thresholds.getLth() * noiseHandle->getNoiseFast(strip, detNoiseRange) + 0.5);
160  theFEDhighThresh =
161  static_cast<int16_t>(thresholds.getHth() * noiseHandle->getNoiseFast(strip, detNoiseRange) + 0.5);
162 
163  adcPrev = -9999;
164  adcNext = -9999;
165  /*
166  If a strip is the last one on the chip
167  set its next neighbor's thresholds to infinity
168  because the FED does not merge clusters across
169  chip boundaries right now
170  */
171  if (strip % 128 == 127) {
172  adcNext = 0;
173  theNextFEDlowThresh = 9999;
174  theNextFEDhighThresh = 9999;
175  } else {
176  adcNext = (in_iter + 1)->adc();
177  SiStripThreshold::Data thresholds_1 = thresholdHandle->getData(strip + 1, detThRange);
178  theNextFEDlowThresh =
179  static_cast<int16_t>(thresholds_1.getLth() * noiseHandle->getNoiseFast(strip + 1, detNoiseRange) + 0.5);
180  theNextFEDhighThresh =
181  static_cast<int16_t>(thresholds_1.getHth() * noiseHandle->getNoiseFast(strip + 1, detNoiseRange) + 0.5);
182  }
183  /*
184  Similarily, for the first strip
185  on a chip
186  */
187  if (strip % 128 == 0) {
188  adcPrev = 0;
189  thePrevFEDlowThresh = 9999;
190  thePrevFEDhighThresh = 9999;
191  } else {
192  adcPrev = (in_iter - 1)->adc();
193  SiStripThreshold::Data thresholds_1 = thresholdHandle->getData(strip - 1, detThRange);
194  thePrevFEDlowThresh =
195  static_cast<int16_t>(thresholds_1.getLth() * noiseHandle->getNoiseFast(strip - 1, detNoiseRange) + 0.5);
196  thePrevFEDhighThresh =
197  static_cast<int16_t>(thresholds_1.getHth() * noiseHandle->getNoiseFast(strip - 1, detNoiseRange) + 0.5);
198  }
199  if (adcNext < adcPrev) {
200  adcMaxNeigh = adcPrev;
201  theNeighFEDlowThresh = thePrevFEDlowThresh;
202  theNeighFEDhighThresh = thePrevFEDhighThresh;
203  } else {
204  adcMaxNeigh = adcNext;
205  theNeighFEDlowThresh = theNextFEDlowThresh;
206  theNeighFEDhighThresh = theNextFEDhighThresh;
207  }
208 
209  //Find adc values for next neighbouring strips
210  adcPrev2 = -9999;
211  adcNext2 = -9999;
212  thePrev2FEDlowThresh = 1;
213  theNext2FEDlowThresh = 1;
214  if (strip % 128 >= 126) {
215  adcNext2 = 0;
216  theNext2FEDlowThresh = 9999;
217  } else if (strip % 128 < 126) {
218  adcNext2 = (in_iter + 2)->adc();
219  theNext2FEDlowThresh = static_cast<int16_t>(thresholdHandle->getData(strip + 2, detThRange).getLth() *
220  noiseHandle->getNoiseFast(strip + 2, detNoiseRange) +
221  0.5);
222  }
223  if (strip % 128 <= 1) {
224  adcPrev2 = 0;
225  thePrev2FEDlowThresh = 9999;
226  } else if (strip % 128 > 1) {
227  adcPrev2 = (in_iter - 2)->adc();
228  thePrev2FEDlowThresh = static_cast<int16_t>(thresholdHandle->getData(strip - 2, detThRange).getLth() *
229  noiseHandle->getNoiseFast(strip - 2, detNoiseRange) +
230  0.5);
231  }
232  //GB 23/6/08: truncation should be done at the very beginning
233  if (isAValidDigi())
234  out.data.push_back(SiStripDigi(strip, truncate(in_iter->adc())));
235  }
236 }
237 
238 void SiStripFedZeroSuppression::fillThresholds_(const uint32_t detID, size_t size) {
239  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detID);
240  SiStripThreshold::Range detThRange = thresholdHandle->getRange(detID);
241 
242  if (highThr_.size() != size) {
243  highThr_.resize(size);
244  lowThr_.resize(size);
245  noises_.resize(size);
246  highThrSN_.resize(size);
247  lowThrSN_.resize(size);
248  }
249 
250  noiseHandle->allNoises(noises_, detNoiseRange);
251  thresholdHandle->allThresholds(lowThrSN_, highThrSN_, detThRange); // thresholds as S/N
252  for (size_t strip = 0; strip < size; ++strip) {
253  float noise = noises_[strip];
254  // uncomment line below to check bluk noise decoding
255  //assert( noise == noiseHandle->getNoiseFast(strip,detNoiseRange) );
256  highThr_[strip] = static_cast<int16_t>(highThrSN_[strip] * noise + 0.5 + 1e-6);
257  lowThr_[strip] = static_cast<int16_t>(lowThrSN_[strip] * noise + 0.5 + 1e-6);
258  // Note: it's a bit wierd, but there are some cases for which 'highThrSN_[strip]*noise' is an exact integer
259  // but due to roundoffs it gets rounded to the integer below if.
260  // Apparently the optimized code inlines differently and this changes the roundoff.
261  // The +1e-6 fixes the problem. [GPetruc]
262  }
263 }
264 
265 void SiStripFedZeroSuppression::suppress(const std::vector<int16_t>& in,
266  uint16_t firstAPV,
268  const uint32_t detID = out.id;
269  size_t size = in.size();
270 #ifdef DEBUG_SiStripZeroSuppression_
271  if (edm::isDebugEnabled())
272  LogTrace("SiStripZeroSuppression")
273  << "[SiStripFedZeroSuppression::suppress] Zero suppression on std::vector<int16_t>: detID " << detID
274  << " size = " << in.size();
275 #endif
276 
277  fillThresholds_(detID, size + firstAPV * 128); // want to decouple this from the other cost
278 
279  std::vector<int16_t>::const_iterator in_iter = in.begin();
280  uint16_t strip = firstAPV * 128;
281  for (; strip < size + firstAPV * 128; ++strip, ++in_iter) {
282  size_t strip_mod_128 = strip & 127;
283 #ifdef DEBUG_SiStripZeroSuppression_
284  if (edm::isDebugEnabled())
285  LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] detID= " << detID
286  << " strip= " << strip << " adc= " << *in_iter;
287 #endif
288  adc = *in_iter;
289 
290  theFEDlowThresh = lowThr_[strip];
291  theFEDhighThresh = highThr_[strip];
292 
293  //Find adc values for neighbouring strips
294 
295  /*
296  If a strip is the last one on the chip
297  set its next neighbor's thresholds to infinity
298  because the FED does not merge clusters across
299  chip boundaries right now
300  */
301 
302  //adcPrev = -9999; // useless, they are set
303  //adcNext = -9999; // in the next lines in any case
304  if (strip_mod_128 == 127) {
305  adcNext = 0;
306  theNextFEDlowThresh = 9999;
307  theNextFEDhighThresh = 9999;
308  } else {
309  adcNext = *(in_iter + 1);
310  theNextFEDlowThresh = lowThr_[strip + 1];
311  theNextFEDhighThresh = highThr_[strip + 1];
312  }
313 
314  /*
315  Similarily, for the first strip
316  on a chip
317  */
318  if (strip_mod_128 == 0) {
319  adcPrev = 0;
320  thePrevFEDlowThresh = 9999;
321  thePrevFEDhighThresh = 9999;
322  } else {
323  adcPrev = *(in_iter - 1);
324  thePrevFEDlowThresh = lowThr_[strip - 1];
325  thePrevFEDhighThresh = highThr_[strip - 1];
326  }
327 
328  if (adcNext < adcPrev) {
329  adcMaxNeigh = adcPrev;
330  theNeighFEDlowThresh = thePrevFEDlowThresh;
331  theNeighFEDhighThresh = thePrevFEDhighThresh;
332  } else {
333  adcMaxNeigh = adcNext;
334  theNeighFEDlowThresh = theNextFEDlowThresh;
335  theNeighFEDhighThresh = theNextFEDhighThresh;
336  }
337 
338  //Find adc values for next neighbouring strips
339  //adcPrev2 = -9999; //
340  //adcNext2 = -9999; // useless to set them here
341  //thePrev2FEDlowThresh = 1; // they are overwritten always in the next 8 lines
342  //theNext2FEDlowThresh = 1; //
343  if (strip_mod_128 >= 126) {
344  adcNext2 = 0;
345  theNext2FEDlowThresh = 9999;
346  //} else if ( strip_mod_128 < 126 ) { // if it's not >= then is <, no need to "if" again
347  } else {
348  adcNext2 = *(in_iter + 2);
349  theNext2FEDlowThresh = lowThr_[strip + 2];
350  }
351  if (strip_mod_128 <= 1) {
352  adcPrev2 = 0;
353  thePrev2FEDlowThresh = 9999;
354  //} else if ( strip_mod_128 > 1 ) { // same as above
355  } else {
356  adcPrev2 = *(in_iter - 2);
357  thePrev2FEDlowThresh = lowThr_[strip - 2];
358  ;
359  }
360 
361  if (isAValidDigi()) {
362 #ifdef DEBUG_SiStripZeroSuppression_
363  if (edm::isDebugEnabled())
364  LogTrace("SiStripZeroSuppression")
365  << "[SiStripFedZeroSuppression::suppress] DetId " << out.id << " strip " << strip << " adc " << *in_iter
366  << " digiCollection size " << out.data.size();
367 #endif
368  //GB 23/6/08: truncation should be done at the very beginning
369  out.push_back(SiStripDigi(strip, (*in_iter < 0 ? 0 : truncate(*in_iter))));
370  }
371  }
372 }
373 
375 #ifdef DEBUG_SiStripZeroSuppression_
376 
377  if (edm::isDebugEnabled()) {
378  LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] "
379  << "\n\t adc " << adc << "\n\t adcPrev " << adcPrev << "\n\t adcNext " << adcNext
380  << "\n\t adcMaxNeigh " << adcMaxNeigh << "\n\t adcPrev2 " << adcPrev2
381  << "\n\t adcNext2 " << adcNext2 << std::endl;
382 
383  LogTrace("SiStripZeroSuppression") << "[SiStripFedZeroSuppression::suppress] "
384  << "\n\t theFEDlowThresh " << theFEDlowThresh << "\n\t theFEDhighThresh "
385  << theFEDhighThresh << "\n\t thePrevFEDlowThresh " << thePrevFEDlowThresh
386  << "\n\t thePrevFEDhighThresh " << thePrevFEDhighThresh
387  << "\n\t theNextFEDlowThresh " << theNextFEDlowThresh
388  << "\n\t theNextFEDhighThresh " << theNextFEDhighThresh
389  << "\n\t theNeighFEDlowThresh " << theNeighFEDlowThresh
390  << "\n\t theNeighFEDhighThresh " << theNeighFEDhighThresh
391  << "\n\t thePrev2FEDlowThresh " << thePrev2FEDlowThresh
392  << "\n\t theNext2FEDlowThresh " << theNext2FEDlowThresh << std::endl;
393  }
394 #endif
395  // Decide if this strip should be accepted.
396  bool accept = false;
397  switch (theFEDalgorithm) {
398  case 1:
399  accept = (adc >= theFEDlowThresh);
400  break;
401  case 2:
402  accept = (adc >= theFEDhighThresh || (adc >= theFEDlowThresh && adcMaxNeigh >= theNeighFEDlowThresh));
403  break;
404  case 3:
405  accept = (adc >= theFEDhighThresh || (adc >= theFEDlowThresh && adcMaxNeigh >= theNeighFEDhighThresh));
406  break;
407  case 4:
408  accept =
409  ((adc >= theFEDhighThresh) //Test for adc>highThresh (same as algorithm 2)
410  || ((adc >= theFEDlowThresh) //Test for adc>lowThresh, with neighbour adc>lowThresh (same as algorithm 2)
411  && (adcMaxNeigh >= theNeighFEDlowThresh)) ||
412  ((adc < theFEDlowThresh) //Test for adc<lowThresh
413  && (((adcPrev >= thePrevFEDhighThresh) //with both neighbours>highThresh
414  && (adcNext >= theNextFEDhighThresh)) ||
415  ((adcPrev >= thePrevFEDhighThresh) //OR with previous neighbour>highThresh and
416  && (adcNext >= theNextFEDlowThresh) //both the next neighbours>lowThresh
417  && (adcNext2 >= theNext2FEDlowThresh)) ||
418  ((adcNext >= theNextFEDhighThresh) //OR with next neighbour>highThresh and
419  && (adcPrev >= thePrevFEDlowThresh) //both the previous neighbours>lowThresh
420  && (adcPrev2 >= thePrev2FEDlowThresh)) ||
421  ((adcNext >= theNextFEDlowThresh) //OR with both next neighbours>lowThresh and
422  && (adcNext2 >= theNext2FEDlowThresh) //both the previous neighbours>lowThresh
423  && (adcPrev >= thePrevFEDlowThresh) && (adcPrev2 >= thePrev2FEDlowThresh)))));
424  break;
425  case 5:
426  accept = adc > 0;
427  break;
428  }
429  return accept;
430 }
mps_fire.i
i
Definition: mps_fire.py:428
SiStripNoises::getNoiseFast
static float getNoiseFast(const uint16_t &strip, const Range &range)
Definition: SiStripNoises.h:68
particleFlowZeroSuppressionECAL_cff.thresholds
thresholds
Definition: particleFlowZeroSuppressionECAL_cff.py:31
edm::DetSet
Definition: DetSet.h:23
SiStripThreshold::Data::getLth
float getLth() const
Definition: SiStripThreshold.h:61
digitizers_cfi.strip
strip
Definition: digitizers_cfi.py:19
SiStripThreshold::Data::getHth
float getHth() const
Definition: SiStripThreshold.h:60
SiStripNoises.h
SiStripNoises::getRange
const Range getRange(const uint32_t detID) const
Definition: SiStripNoises.cc:34
SiStripNoisesRcd
Definition: SiStripCondDataRecords.h:40
SiStripFedZeroSuppression::suppress
void suppress(const std::vector< SiStripDigi > &in, std::vector< SiStripDigi > &selectedSignal, uint32_t detId, edm::ESHandle< SiStripNoises > &, edm::ESHandle< SiStripThreshold > &)
Definition: SiStripFedZeroSuppression.cc:32
suppress
Definition: GenABIO.cc:168
ecalLiteDTU::adc
constexpr int adc(sample_type sample)
get the ADC sample (12 bits)
Definition: EcalLiteDTUSample.h:12
accept
bool accept(const edm::Event &event, const edm::TriggerResults &triggerTable, const std::string &triggerPath)
Definition: TopDQMHelpers.h:30
edm::EventSetup::get
T get() const
Definition: EventSetup.h:80
SiStripFedZeroSuppression.h
SiStripNoises::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripNoises.h:47
SiStripThreshold::getRange
const Range getRange(const uint32_t &detID) const
Definition: SiStripThreshold.cc:38
SiStripFedZeroSuppression::isAValidDigi
bool isAValidDigi()
Definition: SiStripFedZeroSuppression.cc:374
edm::ESHandle< SiStripNoises >
SiStripThreshold::getData
SiStripThreshold::Data getData(const uint16_t &strip, const Range &range) const
Definition: SiStripThreshold.cc:70
SiStripThreshold.h
SiStripThresholdRcd
Definition: SiStripCondDataRecords.h:50
recoMuon::in
Definition: RecoMuonEnumerators.h:6
hgcalDigitizer_cfi.noise
noise
Definition: hgcalDigitizer_cfi.py:155
SiStripThreshold::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripThreshold.h:106
SiStripThresholdRcd.h
edm::EventSetup
Definition: EventSetup.h:57
get
#define get
std
Definition: JetResolutionObject.h:76
SiStripFedZeroSuppression::init
void init(const edm::EventSetup &es)
Definition: SiStripFedZeroSuppression.cc:12
SiStripFedZeroSuppression::fillThresholds_
void fillThresholds_(const uint32_t detID, size_t size)
Definition: SiStripFedZeroSuppression.cc:238
edm::isDebugEnabled
bool isDebugEnabled()
Definition: MessageLogger.cc:12
MillePedeFileConverter_cfg.out
out
Definition: MillePedeFileConverter_cfg.py:31
SiStripDigi
A Digi for the silicon strip detector, containing both strip and adc information, and suitable for st...
Definition: SiStripDigi.h:12
SiStripNoisesRcd.h
LogTrace
#define LogTrace(id)
Definition: MessageLogger.h:224
SiStripThreshold::Data
Definition: SiStripThreshold.h:44
findQualityFiles.size
size
Write out results.
Definition: findQualityFiles.py:443
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
edm::DetSet::const_iterator
collection_type::const_iterator const_iterator
Definition: DetSet.h:31