CMS 3D CMS Logo
SiStripAPVRestorer.cc
Go to the documentation of this file.
1 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripAPVRestorer.h"
2 #include "FWCore/ParameterSet/interface/ParameterSet.h"
3 
4 #include <cmath>
5 #include <iostream>
6 #include <algorithm>
7 
8 
9 SiStripAPVRestorer::SiStripAPVRestorer(const edm::ParameterSet& conf):
10  quality_cache_id(-1), noise_cache_id(-1), pedestal_cache_id(-1),
11  ForceNoRestore_(conf.getParameter<bool>("ForceNoRestore")),
12  SelfSelectRestoreAlgo_(conf.getParameter<bool>("SelfSelectRestoreAlgo")),
13  InspectAlgo_(conf.getParameter<std::string>("APVInspectMode")),
14  RestoreAlgo_(conf.getParameter<std::string>("APVRestoreMode")),
15  useRealMeanCM_(conf.getParameter<bool>("useRealMeanCM")),
16  fraction_(conf.getParameter<double>("Fraction")),
17  deviation_(conf.getParameter<uint32_t>("Deviation")),
18  restoreThreshold_(conf.getParameter<double>("restoreThreshold")),
19  DeltaCMThreshold_(conf.getParameter<uint32_t>("DeltaCMThreshold")),
20  nSigmaNoiseDerTh_(conf.getParameter<uint32_t>("nSigmaNoiseDerTh")),
21  consecThreshold_(conf.getParameter<uint32_t>("consecThreshold")),
22  hitStripThreshold_(conf.getParameter<uint32_t>("hitStripThreshold")),
23  nSmooth_(conf.getParameter<uint32_t>("nSmooth")),
24  minStripsToFit_(conf.getParameter<uint32_t>("minStripsToFit")),
25  distortionThreshold_(conf.getParameter<uint32_t>("distortionThreshold")),
26  CutToAvoidSignal_(conf.getParameter<double>("CutToAvoidSignal")),
27  nSaturatedStrip_(conf.getParameter<uint32_t>("nSaturatedStrip")),
28  ApplyBaselineCleaner_(conf.getParameter<bool>("ApplyBaselineCleaner")),
29  slopeX_(conf.getParameter<int32_t>("slopeX")),
30  slopeY_(conf.getParameter<int32_t>("slopeY")),
31  CleaningSequence_(conf.getParameter<uint32_t>("CleaningSequence")),
32  ApplyBaselineRejection_(conf.getParameter<bool>("ApplyBaselineRejection")),
33  MeanCM_(conf.getParameter<int32_t>("MeanCM")),
34  filteredBaselineMax_(conf.getParameter<double>("filteredBaselineMax")),
35  filteredBaselineDerivativeSumSquare_(conf.getParameter<double>("filteredBaselineDerivativeSumSquare"))
36 
37 {
38  apvFlags_.clear();
39  median_.clear();
40  SmoothedMaps_.clear();
41  BaselineMap_.erase(BaselineMap_.begin(), BaselineMap_.end());
42 }
43 
44 
45 void SiStripAPVRestorer::init(const edm::EventSetup& es){
46  uint32_t n_cache_id = es.get<SiStripNoisesRcd>().cacheIdentifier();
47  uint32_t q_cache_id = es.get<SiStripQualityRcd>().cacheIdentifier();
48  uint32_t p_cache_id = es.get<SiStripPedestalsRcd>().cacheIdentifier();
49 
50  if(n_cache_id != noise_cache_id) {
51  es.get<SiStripNoisesRcd>().get( noiseHandle );
52  noise_cache_id = n_cache_id;
53  } else {
54  noise_cache_id = n_cache_id;
55  }
56  if(q_cache_id != quality_cache_id) {
57  es.get<SiStripQualityRcd>().get( qualityHandle );
58  quality_cache_id = q_cache_id;
59  }else {
60  quality_cache_id = q_cache_id;
61  }
62 
63  if(p_cache_id != pedestal_cache_id) {
64  es.get<SiStripPedestalsRcd>().get( pedestalHandle );
65  pedestal_cache_id = p_cache_id;
66  }else {
67  pedestal_cache_id = p_cache_id;
68  }
69 
70 }
71 
72 
73 int16_t SiStripAPVRestorer::InspectAndRestore( const uint32_t& detId, const uint16_t& firstAPV, std::vector<int16_t>& rawDigisPedSubtracted, std::vector<int16_t>& processedRawDigi, const std::vector< std::pair<short,float> >& vmedians ){
74  int16_t nAPVFlagged = this->inspect(detId, firstAPV, rawDigisPedSubtracted, vmedians);
75  this->restore(firstAPV, processedRawDigi);
76  return nAPVFlagged;
77 }
78 
79 
80 int16_t SiStripAPVRestorer::inspect( const uint32_t& detId, const uint16_t& firstAPV, std::vector<int16_t>& digis, const std::vector< std::pair<short,float> >& vmedians) {
81 
82  detId_ = detId;
83 
84  apvFlagsBool_.clear();
85  apvFlagsBoolOverride_.clear();
86  apvFlagsBoolOverride_.insert(apvFlagsBoolOverride_.begin(), 6, false);
87  apvFlags_.clear();
88  apvFlags_.insert(apvFlags_.begin(), 6, "");
89  median_.clear();
90  median_.insert(median_.begin(), 6, -999);
91  badAPVs_.clear();
92  badAPVs_.insert(badAPVs_.begin(), 6, false);
93  SmoothedMaps_.erase(SmoothedMaps_.begin(), SmoothedMaps_.end());
94  BaselineMap_.erase(BaselineMap_.begin(), BaselineMap_.end());
95 
96  for(size_t i=0; i< vmedians.size(); ++i){
97  short APV = vmedians[i].first;
98  median_[APV]= vmedians[i].second;
99  badAPVs_[APV] = qualityHandle->IsApvBad(detId_, APV);
100  }
101 
102  if(InspectAlgo_=="BaselineFollower") return this->BaselineFollowerInspect(firstAPV, digis);
103  if(InspectAlgo_=="AbnormalBaseline") return this->AbnormalBaselineInspect(firstAPV, digis);
104  if(InspectAlgo_=="Null") return this->NullInspect(firstAPV, digis);
105  if(InspectAlgo_=="BaselineAndSaturation") return this->BaselineAndSaturationInspect(firstAPV, digis);
106  throw cms::Exception("Unregistered Inspect Algorithm") << "SiStripAPVRestorer possibilities: (Null), (AbnormalBaseline),(BaselineFollower)";
107 
108 }
109 
110 
111 void SiStripAPVRestorer::restore(const uint16_t& firstAPV, std::vector<int16_t>& digis ) {
112 
113  if(ForceNoRestore_) return;
114 
115  for( uint16_t APV=firstAPV; APV< digis.size()/128 + firstAPV; ++APV){
116  std::string algoToUse = *( apvFlags_.begin() + APV );
117 
118  if ( algoToUse != ""){
119  if(!SelfSelectRestoreAlgo_) algoToUse = RestoreAlgo_;
120 
121  if(algoToUse=="Flat"){
122  this->FlatRestore(APV, firstAPV, digis);
123  }else if(algoToUse=="BaselineFollower"){
124  this->BaselineFollowerRestore(APV, firstAPV, median_[APV], digis);
125  }else{
126  throw cms::Exception("Unregistered Restore Algorithm") << "SiStripAPVRestorer possibilities: (Flat), (BaselineFollower)";
127  }
128 
129 
130  }
131  }
132 
133 }
134 
135 
136 //Inspect method implementation ========================================================================================================================================================================
137 //======================================================================================================================================================================================================
138 //======================================================================================================================================================================================================
139 
140 template<typename T>
141 inline
142 int16_t SiStripAPVRestorer::BaselineFollowerInspect(const uint16_t& firstAPV, std::vector<T>& digis){
143  std::vector<T> singleAPVdigi;
144  int16_t nAPVflagged = 0;
145 
146  CMMap::iterator itCMMap;
147  if(useRealMeanCM_) itCMMap = MeanCMmap_.find(detId_);
148 
149  for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
150 
151  DigiMap smoothedmap;
152  smoothedmap.erase(smoothedmap.begin(), smoothedmap.end());
153 
154  if(!badAPVs_[APV]){
155  float MeanAPVCM = MeanCM_;
156  if(useRealMeanCM_&&itCMMap!= MeanCMmap_.end()) MeanAPVCM =(itCMMap->second)[APV];
157 
158  singleAPVdigi.clear();
159  for(int16_t strip = (APV-firstAPV)*128; strip < (APV-firstAPV+1)*128; ++strip){
160  singleAPVdigi.push_back(digis[strip]);
161  }
162 
163 
164  float DeltaCM = median_[APV] - MeanAPVCM;
165 
166  //std::cout << "Delta CM: " << DeltaCM << " CM: " << median_[APV] << " detId " << (uint32_t) detId_ << std::endl;
167  if(DeltaCM < 0 && std::abs(DeltaCM) > DeltaCMThreshold_){
168 
169  bool isFlat = FlatRegionsFinder(singleAPVdigi,smoothedmap,APV);
170  if(!isFlat){
171  apvFlags_[APV]= "BaselineFollower"; //specify any algo to make the restore
172  nAPVflagged++;
173  }
174  }
175 
176  }
177  SmoothedMaps_.insert(SmoothedMaps_.end(), std::pair<uint16_t, DigiMap>(APV, smoothedmap));
178  }
179 
180  return nAPVflagged;
181 }
182 
183 //======================================================================================================================================================================================================
184 template<typename T>
185 inline
186 int16_t SiStripAPVRestorer::BaselineAndSaturationInspect(const uint16_t& firstAPV, std::vector<T>& digis){
187  std::vector<T> singleAPVdigi;
188  singleAPVdigi.clear();
189 
190 
191  int16_t nAPVflagged = 0;
192 
193  CMMap::iterator itCMMap;
194  if(useRealMeanCM_) itCMMap = MeanCMmap_.find(detId_);
195 
196  for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
197  apvFlags_.push_back( "" );
198  if(!badAPVs_[APV]){
199  float MeanAPVCM = MeanCM_;
200  if(useRealMeanCM_&&itCMMap!= MeanCMmap_.end()) MeanAPVCM =(itCMMap->second)[APV];
201 
202  singleAPVdigi.clear();
203 
204  uint16_t nSatStrip =0;
205  for(int16_t strip = (APV-firstAPV)*128; strip < (APV-firstAPV+1)*128; ++strip){
206  singleAPVdigi.push_back(digis[strip]);
207  if(digis[strip] >=1023) ++nSatStrip;
208  }
209 
210  float DeltaCM = median_[APV] -MeanAPVCM;
211 
212 
213  if(DeltaCM < 0 && std::abs(DeltaCM) > DeltaCMThreshold_&&nSatStrip>= nSaturatedStrip_){
214  apvFlags_[APV] = RestoreAlgo_; //specify any algo to make the restore
215  nAPVflagged++;
216  }
217  }
218  }
219 
220  return nAPVflagged;
221 }
222 
223 //======================================================================================================================================================================================================
224 template<typename T>
225 inline
226 int16_t SiStripAPVRestorer::AbnormalBaselineInspect( const uint16_t& firstAPV, std::vector<T>& digis){
227 
228  SiStripQuality::Range detQualityRange = qualityHandle->getRange(detId_);
229 
230  typename std::vector<T>::iterator fs;
231 
232  int16_t nAPVflagged=0;
233 
234  CMMap::iterator itCMMap;
235  if(useRealMeanCM_) itCMMap = MeanCMmap_.find(detId_);
236 
237 
238  int devCount = 0, qualityCount = 0, minstrip = 0;
239  for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
240  apvFlags_.push_back( "" );
241  if(!badAPVs_[APV]){
242  float MeanAPVCM = MeanCM_;
243  if(useRealMeanCM_&&itCMMap!= MeanCMmap_.end()) MeanAPVCM =(itCMMap->second)[APV];
244  for (uint16_t istrip=APV*128; istrip<(APV+1)*128; ++istrip){
245  fs = digis.begin() + istrip-firstAPV*128;
246  if ( !qualityHandle->IsStripBad(detQualityRange,istrip) ){
247  qualityCount++;
248  if ( std::abs((int) *fs - MeanAPVCM) > (int)deviation_ ){
249  devCount++;
250  minstrip = std::min((int) *fs, minstrip);
251  }
252  }
253  }
254 
255  if( devCount > fraction_ * qualityCount ) {
256  apvFlags_[APV] = RestoreAlgo_; //specify any algo to make the restore
257  nAPVflagged++;
258  }
259  }
260  }
261 
262  return nAPVflagged;
263 
264 }
265 
266 
267 //======================================================================================================================================================================================================
268 template<typename T>
269 inline
270 int16_t SiStripAPVRestorer::NullInspect(const uint16_t& firstAPV, std::vector<T>& digis){
271 
272  SiStripQuality::Range detQualityRange = qualityHandle->getRange(detId_);
273 
274  typename std::vector<T>::iterator fs;
275 
276  int16_t nAPVflagged = 0;
277 
278  for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
279  apvFlags_.push_back( "" );
280  if(!badAPVs_[APV]){
281  int zeroCount = 0, qualityCount = 0;
282  for (uint16_t istrip=APV*128; istrip<(APV+1)*128; ++istrip){
283  fs = digis.begin() + istrip-firstAPV*128;
284  if ( !qualityHandle->IsStripBad(detQualityRange,istrip) ){
285  qualityCount++;
286  if ( (int) *fs < 1 ) zeroCount++;
287  }
288  }
289 
290  if( zeroCount > restoreThreshold_ * qualityCount ) {
291  apvFlags_[APV] = RestoreAlgo_; //specify any algo to make the restore
292  nAPVflagged++;
293  }
294  }
295  }
296 
297  return nAPVflagged;
298 
299 }
300 
301 
302 
303 //Restore method implementation ========================================================================================================================================================================
304 //======================================================================================================================================================================================================
305 //======================================================================================================================================================================================================
306 
307 inline
308 void SiStripAPVRestorer::BaselineFollowerRestore(const uint16_t& APVn, const uint16_t& firstAPV, const float& median, std::vector<int16_t>& digis){
309  //typename std::vector<T>::iterator firstStrip(digis.begin() + APVn*128), lastStrip(firstStrip + 128), actualStrip;
310 
311 
312  std::vector<int16_t> baseline;
313  baseline.clear();
314  baseline.insert(baseline.begin(),128, 0);
315 
316 
317 
318 
319  //============================= Find Flat Regions & Interpolating the baseline & subtracting the baseline =================
320 
321  if(SmoothedMaps_.size()){
322  std::map<uint16_t, DigiMap >::iterator itSmootedMap = SmoothedMaps_.find(APVn);
323  this->BaselineFollower(itSmootedMap->second, baseline, median);
324  } else {
325  //median=0;
326  DigiMap smoothedpoints;
327  std::vector<int16_t> singleAPVdigi;
328  singleAPVdigi.clear();
329  for(int16_t strip = (APVn-firstAPV)*128; strip < (APVn-firstAPV+1)*128; ++strip) singleAPVdigi.push_back(digis[strip]);
330  this->FlatRegionsFinder(singleAPVdigi,smoothedpoints, APVn);
331  this->BaselineFollower(smoothedpoints, baseline, median);
332 
333  }
334 
335  if(ApplyBaselineRejection_){
336  if(CheckBaseline(baseline)) apvFlagsBoolOverride_[APVn] = true;
337  }
338 
339  //============================= subtracting the baseline =============================================
340 
341  for(int16_t itStrip= 0 ; itStrip< 128; ++itStrip){
342  digis[(APVn-firstAPV)*128+itStrip] -= baseline[itStrip] - median;
343  }
344 
345 
346  //============================= storing baseline to the map =============================================
347  BaselineMap_.insert(BaselineMap_.end(), std::pair< uint16_t, std::vector < int16_t> >(APVn, baseline));
348 
349 }
350 
351 
352 //======================================================================================================================================================================================================
353 inline
354 void SiStripAPVRestorer::FlatRestore(const uint16_t& APVn, const uint16_t& firstAPV, std::vector<int16_t>& digis ){
355 
356  std::vector<int16_t> baseline;
357  baseline.clear();
358  baseline.insert(baseline.begin(),128, 150);
359  baseline[0]=0; baseline[127]=0;
360  BaselineMap_.insert(BaselineMap_.end(), std::pair< uint16_t, std::vector < int16_t> >(APVn, baseline));
361 
362  for(int16_t itStrip= 0 ; itStrip< 128; ++itStrip){
363  digis[(APVn-firstAPV)*128+itStrip] = baseline[itStrip];
364  }
365 
366 
367 }
368 
369 
370 
371 //Baseline calculation implementation ==================================================================================================================================================================
372 //======================================================================================================================================================================================================
373 //======================================================================================================================================================================================================
374 
375 bool inline SiStripAPVRestorer::FlatRegionsFinder(const std::vector<int16_t>& adcs, DigiMap& smoothedpoints, const uint16_t& APVn){
376  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detId_);
377 
378  DigiMap consecpoints;
379  DigiMapIter itConsecpoints, itSmoothedpoints;
380  consecpoints.erase(consecpoints.begin(), consecpoints.end());
381  smoothedpoints.erase(smoothedpoints.begin(), smoothedpoints.end());
382 
383 
384  //============================= Height above local minimum ===============================
385  std::vector<float> adcsLocalMinSubtracted;
386  adcsLocalMinSubtracted.clear();
387  adcsLocalMinSubtracted.insert(adcsLocalMinSubtracted.begin(), 128,0);
388  for(uint32_t istrip=0; istrip<128; ++istrip) {
389  float localmin = 999.9;
390  for(uint16_t jstrip=std::max(0,(int)(istrip-nSmooth_/2)); jstrip<std::min(128,(int)(istrip+nSmooth_/2)); ++jstrip) {
391  float nextvalue = adcs[jstrip];
392  if(nextvalue < localmin) localmin=nextvalue;
393  }
394  adcsLocalMinSubtracted[istrip] = adcs[istrip] - localmin;
395  }
396 
397 
398  //============================= Find regions with stable slopes ========================
399  std::vector<uint16_t> nConsStrip;
400  nConsStrip.clear();
401 
402  //Creating maps with all the neighborhood strip and putting in a nCosntStip vector how many we have
403  uint16_t consecStrips=0;
404  for(uint32_t istrip=0; istrip<128; ++istrip) {
405  int16_t adc = adcs[istrip];
406 
407  //if( adcsLocalMinSubtracted[istrip] < nSigmaNoiseDerTh_ * (float)noiseHandle->getNoise(istrip+APVn*128,detNoiseRange) && (adc - median) < hitStripThreshold_){
408  if( adcsLocalMinSubtracted[istrip] < nSigmaNoiseDerTh_ * (float)noiseHandle->getNoiseFast(istrip+APVn*128,detNoiseRange)){
409  consecpoints.insert(consecpoints.end(), std::pair<uint16_t, int16_t >(istrip, adc));
410  ++consecStrips;
411  }else if (consecStrips >0){
412  nConsStrip.push_back(consecStrips);
413  consecStrips = 0;
414  }
415  }
416 
417  //to cope with the last flat region of the APV
418  if(consecStrips >0) nConsStrip.push_back(consecStrips);
419 
420  //removing from the map the fist and last points in wide flat regions and erasing from the map too small regions
421  itConsecpoints = consecpoints.begin();
422  float MinSmoothValue=20000., MaxSmoothValue=0.;
423  for(std::vector<uint16_t>::iterator itnConsStrip = nConsStrip.begin(); itnConsStrip < nConsStrip.end(); ++itnConsStrip){
424 
425  consecStrips = *itnConsStrip;
426  if(consecStrips >=consecThreshold_){
427  ++itConsecpoints; //skipping first point
428  uint16_t nFirstStrip = itConsecpoints->first;
429  uint16_t nLastStrip;
430  float smoothValue = 0.0;
431  float stripCount =1;
432  for(uint16_t n =0; n < consecStrips-2; ++n){
433  smoothValue += itConsecpoints->second;
434  if(stripCount == consecThreshold_){
435  smoothValue /= (float)stripCount;
436  nLastStrip = nFirstStrip + stripCount -1;
437  smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nFirstStrip, smoothValue));
438  smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nLastStrip, smoothValue));
439  if(smoothValue > MaxSmoothValue) MaxSmoothValue = smoothValue;
440  if(smoothValue < MinSmoothValue) MinSmoothValue = smoothValue;
441  nFirstStrip = nLastStrip+1;
442  smoothValue=0;
443  stripCount=0;
444  }
445  ++stripCount;
446  ++itConsecpoints;
447  }
448  ++itConsecpoints; //and putting the pointer to the new seies of point
449 
450  if(stripCount>1) {
451  //if(smoothValue>0){
452  --stripCount;
453  smoothValue /= (float)(stripCount);
454  nLastStrip = nFirstStrip + stripCount -1;
455  smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nFirstStrip, smoothValue));
456  smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nLastStrip, smoothValue));
457  if(smoothValue > MaxSmoothValue) MaxSmoothValue = smoothValue;
458  if(smoothValue < MinSmoothValue) MinSmoothValue = smoothValue;
459  }
460  } else{
461  for(int n =0; n< consecStrips ; ++n) ++itConsecpoints;
462  }
463  }
464 
465 
466  if( (MaxSmoothValue-MinSmoothValue) > distortionThreshold_){
467  if(ApplyBaselineCleaner_) this->BaselineCleaner(adcs, smoothedpoints, APVn);
468  return false;
469  }
470  return true;
471 }
472 
473 
474 void inline SiStripAPVRestorer::BaselineCleaner(const std::vector<int16_t>& adcs, DigiMap& smoothedpoints, const uint16_t& APVn){
475 
476  if(CleaningSequence_==0) { //default sequence used up to now
477  this->Cleaner_HighSlopeChecker(smoothedpoints);
478  this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
479  }else if(CleaningSequence_==1){
480  this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
481  this->Cleaner_HighSlopeChecker(smoothedpoints);
482  this->Cleaner_MonotonyChecker(smoothedpoints);
483  }else if(CleaningSequence_==2){
484  this->Cleaner_HighSlopeChecker(smoothedpoints);
485  }else if(CleaningSequence_==3){
486  this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
487  this->Cleaner_HighSlopeChecker(smoothedpoints);
488  }else{
489  this->Cleaner_HighSlopeChecker(smoothedpoints);
490  this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
491  }
492 
493 }
494 
495 
496 void inline SiStripAPVRestorer::Cleaner_MonotonyChecker(DigiMap& smoothedpoints){
497 //Removing points without monotony
498  //--------------------------------------------------------------------------------------------------
499  if(smoothedpoints.size() < 3) return;
500  DigiMapIter itSmoothedpoints, itSmoothedpointsNext, itSmoothedpointsNextNext, itSmoothedpointsBegin, itSmoothedpointsEnd;
501 
502  itSmoothedpoints=smoothedpoints.begin();
503  while (smoothedpoints.size() > 3 && itSmoothedpoints != --(--(smoothedpoints.end()))) { //while we are not at the last point
504  // get info about current and next points
505  itSmoothedpointsNext = itSmoothedpoints;
506  ++itSmoothedpointsNext;
507  itSmoothedpointsNextNext = itSmoothedpointsNext;
508  ++itSmoothedpointsNextNext;
509  float adc1 = itSmoothedpoints->second;
510  float adc2 = itSmoothedpointsNext->second;
511  float adc3 = itSmoothedpointsNextNext->second;
512 
513  if((adc2-adc1) > hitStripThreshold_ && (adc2-adc3) > hitStripThreshold_){
514  smoothedpoints.erase(itSmoothedpointsNext);
515  }else {
516  ++itSmoothedpoints;
517  }
518 
519  }
520 }
521 
522 void inline SiStripAPVRestorer::Cleaner_LocalMinimumAdder(const std::vector<int16_t>& adcs, DigiMap& smoothedpoints, const uint16_t& APVn){
523  SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detId_);
524  //inserting extra point is case of local minimum
525  //--------------------------------------------------------------------------------------------------
526  // these should be reset now for the point-insertion that follows
527 
528  DigiMapIter itSmoothedpoints, itSmoothedpointsNext, itSmoothedpointsBegin, itSmoothedpointsEnd;
529  if(smoothedpoints.size() >= 2){
530 
531  itSmoothedpointsBegin = smoothedpoints.begin();
532  itSmoothedpointsEnd = --(smoothedpoints.end());
533  for(itSmoothedpoints = itSmoothedpointsBegin; itSmoothedpoints != itSmoothedpointsEnd; ++itSmoothedpoints){
534  itSmoothedpointsNext = itSmoothedpoints;
535  ++itSmoothedpointsNext;
536  float strip1 = itSmoothedpoints->first;
537  float strip2 = itSmoothedpointsNext->first;
538  float adc1 = itSmoothedpoints->second;
539  float adc2 = itSmoothedpointsNext->second;
540  float m = (adc2 -adc1)/(strip2 -strip1);
541 
542  //2,4
543  if((strip2 - strip1) >slopeX_ && std::abs(adc1 -adc2) >slopeY_){
544  float itStrip = 1;
545  float strip = itStrip + strip1;
546  while(strip < strip2){
547 
548  float adc = adcs[strip];
549  if( adc < (adc1 + m * itStrip - 2 * (float)noiseHandle->getNoiseFast(strip+APVn*128,detNoiseRange))){
550  //std::cout << "applying correction strip: " << strip + APVn*128 << " adc " << adc << " detId: " << detId_ << std::endl;
551  smoothedpoints.insert(itSmoothedpointsNext, std::pair<uint16_t, int16_t >(strip,adc));
552  ++itSmoothedpoints;
553  ++itSmoothedpointsNext;
554  itSmoothedpointsEnd = --(smoothedpoints.end());
555  }
556  ++itStrip;
557  ++strip;
558  }
559 
560 
561  }
562  }
563  }
564 
565 
566  itSmoothedpointsBegin = smoothedpoints.begin();
567  itSmoothedpointsEnd = --(smoothedpoints.end());
568  uint16_t firstStripFlat = itSmoothedpointsBegin->first;
569  uint16_t lastStripFlat = itSmoothedpointsEnd->first;
570  int16_t firstStripFlatADC= itSmoothedpointsBegin->second;
571  int16_t lastStripFlatADC= itSmoothedpointsEnd->second;
572 
573  itSmoothedpoints = itSmoothedpointsBegin;
574  if(firstStripFlat >3){
575  float strip = 0;
576  while(strip < firstStripFlat){
577  float adc = adcs[strip];
578  if( adc < ( firstStripFlatADC - 2 * (float)noiseHandle->getNoiseFast(strip+APVn*128,detNoiseRange))){
579  smoothedpoints.insert(itSmoothedpoints, std::pair<uint16_t, int16_t >(strip,adc));
580  ++itSmoothedpoints;
581  }
582  ++strip;
583  }
584  }
585 
586  itSmoothedpoints = itSmoothedpointsEnd;
587  if(lastStripFlat <125){
588  float strip = lastStripFlat+1;
589  while(strip < 128){
590  float adc = adcs[strip];
591  if( adc < ( lastStripFlatADC - 2 * (float)noiseHandle->getNoiseFast(strip+APVn*128,detNoiseRange))){
592  smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(strip,adc));
593  }
594  ++strip;
595  }
596  }
597 }
598 
599 
600 void inline SiStripAPVRestorer::Cleaner_HighSlopeChecker(DigiMap& smoothedpoints){
601  //Removing points in the slope is too high
602  //--------------------------------------------------------------------------------------------------
603 
604  if(smoothedpoints.size() < 4) return;
605  DigiMapIter itSmoothedpoints, itSmoothedpointsNext, itSmoothedpointsBegin, itSmoothedpointsEnd;
606  itSmoothedpoints=smoothedpoints.begin();
607  while (smoothedpoints.size() >2 && itSmoothedpoints != --(smoothedpoints.end()) ) { //while we are not at the last point
608  //if(smoothedpoints.size() <2) break;
609  // get info about current and next points
610  itSmoothedpointsNext = itSmoothedpoints;
611  ++itSmoothedpointsNext;
612  float strip1 = itSmoothedpoints->first;
613  float strip2 = itSmoothedpointsNext->first;
614  float adc1 = itSmoothedpoints->second;
615  float adc2 = itSmoothedpointsNext->second;
616  float m = (adc2 -adc1)/(strip2 -strip1);
617 
618  if (m>2) { // in case of large positive slope, remove next point and try again from same current point
619  smoothedpoints.erase(itSmoothedpointsNext);
620  } else if (m<-2) { // in case of large negative slope, remove current point and either...
621  // move to next point if we have reached the beginning (post-increment to avoid invalidating pointer during erase) or...
622  if(itSmoothedpoints==smoothedpoints.begin()) smoothedpoints.erase(itSmoothedpoints++);
623  // try again from the previous point if we have not reached the beginning
624  else smoothedpoints.erase(itSmoothedpoints--);
625  } else { // in case of a flat enough slope, continue on to the next point
626  itSmoothedpoints++;
627  }
628 
629  }
630 }
631 
632 void inline SiStripAPVRestorer::BaselineFollower(DigiMap& smoothedpoints, std::vector<int16_t>& baseline, const float& median){
633 
634  baseline.clear();
635  DigiMapIter itSmoothedpoints;
636 
637 
638  //if not enough points
639  if(smoothedpoints.size() < minStripsToFit_){
640  baseline.insert(baseline.begin(),128, median);
641  } else {
642  baseline.insert(baseline.begin(),128, 0);
643 
644  DigiMapIter itSmoothedpointsBegin, itSmoothedpointsEnd;
645  itSmoothedpointsBegin = smoothedpoints.begin();
646  itSmoothedpointsEnd = --(smoothedpoints.end());
647 
648 
649  uint16_t firstStripFlat = itSmoothedpointsBegin->first;
650  uint16_t lastStripFlat = itSmoothedpointsEnd->first;
651  int16_t firstStripFlatADC= itSmoothedpointsBegin->second;
652  int16_t lastStripFlatADC= itSmoothedpointsEnd->second;
653 
654  //adding here the costant line at the extremities
655  baseline.erase(baseline.begin(), baseline.begin()+firstStripFlat);
656  baseline.insert(baseline.begin(), firstStripFlat, firstStripFlatADC);
657 
658  baseline.erase(baseline.begin()+lastStripFlat, baseline.end());
659  baseline.insert(baseline.end(), 128 - lastStripFlat, lastStripFlatADC);
660 
661 
662  //IMPORTANT: the itSmoothedpointsEnd should be at least smaller than smoothedpoints.end() -1
663  for(itSmoothedpoints = itSmoothedpointsBegin; itSmoothedpoints != itSmoothedpointsEnd; ++itSmoothedpoints){
664  DigiMapIter itSmoothedpointsNext = itSmoothedpoints;
665  ++itSmoothedpointsNext;
666  float strip1 = itSmoothedpoints->first;
667  float strip2 = itSmoothedpointsNext->first;
668  float adc1 = itSmoothedpoints->second;
669  float adc2 = itSmoothedpointsNext->second;
670 
671  baseline[strip1] = adc1;
672  baseline[strip2] = adc2;
673  float m = (adc2 -adc1)/(strip2 -strip1);
674  uint16_t itStrip = strip1 +1;
675  float stripadc = adc1 + m;
676  while(itStrip < strip2){
677  baseline[itStrip] = stripadc;
678  ++itStrip;
679  stripadc+=m;
680  }
681 
682  }
683 
684  }
685 }
686 
687 
688 bool SiStripAPVRestorer::CheckBaseline(const std::vector<int16_t> &baseline) const
689 {
690  // The Savitzky-Golay (S-G) filter of any left length nL, right
691  // length nR, and order m, and with an optional opt equals the
692  // derivative order (0 for the magnitude, 1 for the first
693  // derivative, etc.) can be calculated using the following
694  // Mathematica code:
695  //
696  // SavitzkyGolay[m_?IntegerQ, {nL_?IntegerQ, nR_?IntegerQ},
697  // opt___] := Module[
698  // {a, d},
699  // d = If[opt === Null, 0, If[IntegerQ[opt], opt, 0]];
700  // a = Table[
701  // If[i == 0 && j == 0, 1, i^j], {i, -nL, nR}, {j, 0,
702  // m}]; (Inverse[Transpose[a].a].Transpose[a])[[d + 1]]];
703  //
704  // The following coefficients can be then calculated by:
705  //
706  // N[Join[Table[SavitzkyGolay[2, {k, 16}], {k, 0, 16}],
707  // Table[SavitzkyGolay[2, {16, k}], {k, 15, 0, -1}]]]
708 
709  // nLR = max(nL, nR)
710  static const size_t savitzky_golay_n_l_r = 16;
711  static const float savitzky_golay_coefficient
712  [2 * savitzky_golay_n_l_r + 1][2 * savitzky_golay_n_l_r + 1] = {
713  { 0.422085, 0.325077, 0.23839, 0.162023, 0.0959752, 0.0402477,
714  -0.00515996, -0.0402477, -0.0650155, -0.0794634, -0.0835913,
715  -0.0773994, -0.0608875, -0.0340557, 0.00309598, 0.0505676,
716  0.108359 },
717  { 0.315789, 0.254902, 0.19969, 0.150155, 0.106295, 0.0681115,
718  0.0356037, 0.00877193, -0.0123839, -0.0278638, -0.0376677,
719  -0.0417957, -0.0402477, -0.0330237, -0.0201238, -0.00154799,
720  0.0227038, 0.0526316 },
721  { 0.234586, 0.198496, 0.165207, 0.134719, 0.107032, 0.0821465,
722  0.0600619, 0.0407784, 0.024296, 0.0106148, -0.000265369,
723  -0.00834439, -0.0136223, -0.0160991, -0.0157747, -0.0126493,
724  -0.00672269, 0.00200501, 0.0135338 },
725  { 0.172078, 0.153076, 0.135099, 0.118148, 0.102221, 0.0873206,
726  0.073445, 0.0605947, 0.0487697, 0.0379699, 0.0281955,
727  0.0194463, 0.0117225, 0.00502392, -0.000649351, -0.00529733,
728  -0.00892003, -0.0115174, -0.0130895, -0.0136364 },
729  { 0.123659, 0.116431, 0.109144, 0.101798, 0.0943921, 0.0869268,
730  0.0794021, 0.0718179, 0.0641743, 0.0564712, 0.0487087,
731  0.0408868, 0.0330054, 0.0250646, 0.0170644, 0.00900473,
732  0.000885613, -0.00729294, -0.0155309, -0.0238283,
733  -0.0321852 },
734  { 0.0859684, 0.0868154, 0.0869565, 0.0863919, 0.0851214,
735  0.0831451, 0.080463, 0.0770751, 0.0729814, 0.0681818,
736  0.0626765, 0.0564653, 0.0495483, 0.0419255, 0.0335968,
737  0.0245624, 0.0148221, 0.00437606, -0.00677583, -0.0186335,
738  -0.0311971, -0.0444664 },
739  { 0.0565217, 0.0628458, 0.0680971, 0.0722756, 0.0753811,
740  0.0774139, 0.0783738, 0.0782609, 0.0770751, 0.0748165,
741  0.071485, 0.0670807, 0.0616036, 0.0550536, 0.0474308,
742  0.0387352, 0.0289667, 0.0181254, 0.00621118, -0.00677583,
743  -0.0208357, -0.0359684, -0.0521739 },
744  { 0.0334615, 0.0434281, 0.0521329, 0.0595759, 0.0657571,
745  0.0706765, 0.0743341, 0.07673, 0.0778641, 0.0777364,
746  0.0763469, 0.0736957, 0.0697826, 0.0646078, 0.0581712,
747  0.0504728, 0.0415126, 0.0312907, 0.0198069, 0.00706142,
748  -0.00694588, -0.022215, -0.0387458, -0.0565385 },
749  { 0.0153846, 0.0276923, 0.0386622, 0.0482943, 0.0565886,
750  0.0635452, 0.0691639, 0.0734448, 0.076388, 0.0779933,
751  0.0782609, 0.0771906, 0.0747826, 0.0710368, 0.0659532,
752  0.0595318, 0.0517726, 0.0426756, 0.0322408, 0.0204682,
753  0.00735786, -0.0070903, -0.0228763, -0.04, -0.0584615 },
754  { 0.001221, 0.0149451, 0.027326, 0.0383639, 0.0480586,
755  0.0564103, 0.0634188, 0.0690842, 0.0734066, 0.0763858,
756  0.078022, 0.078315, 0.077265, 0.0748718, 0.0711355,
757  0.0660562, 0.0596337, 0.0518681, 0.0427595, 0.0323077,
758  0.0205128, 0.00737485, -0.00710623, -0.0229304, -0.0400977,
759  -0.0586081 },
760  { -0.00985222, 0.00463138, 0.0178098, 0.029683, 0.0402509,
761  0.0495137, 0.0574713, 0.0641236, 0.0694708, 0.0735127,
762  0.0762494, 0.0776809, 0.0778073, 0.0766284, 0.0741442,
763  0.0703549, 0.0652604, 0.0588607, 0.0511557, 0.0421456,
764  0.0318302, 0.0202097, 0.0072839, -0.00694708, -0.0224833,
765  -0.0393247, -0.0574713 },
766  { -0.0184729, -0.00369458, 0.00984169, 0.0221359, 0.0331881,
767  0.0429982, 0.0515662, 0.0588923, 0.0649762, 0.0698181,
768  0.073418, 0.0757758, 0.0768915, 0.0767652, 0.0753968,
769  0.0727864, 0.0689339, 0.0638394, 0.0575028, 0.0499242,
770  0.0411035, 0.0310408, 0.019736, 0.00718917, -0.00659972,
771  -0.0216307, -0.0379037, -0.0554187 },
772  { -0.025139, -0.0103925, 0.00318873, 0.0156046, 0.0268552,
773  0.0369405, 0.0458605, 0.0536151, 0.0602045, 0.0656285,
774  0.0698872, 0.0729806, 0.0749086, 0.0756714, 0.0752688,
775  0.0737009, 0.0709677, 0.0670692, 0.0620054, 0.0557763,
776  0.0483818, 0.039822, 0.0300969, 0.0192065, 0.0071508,
777  -0.00607024, -0.0204566, -0.0360083, -0.0527253 },
778  { -0.0302419, -0.0157536, -0.00234785, 0.00997537, 0.021216,
779  0.0313741, 0.0404497, 0.0484427, 0.0553532, 0.0611811,
780  0.0659264, 0.0695892, 0.0721695, 0.0736672, 0.0740823,
781  0.0734149, 0.0716649, 0.0688324, 0.0649174, 0.0599198,
782  0.0538396, 0.0466769, 0.0384316, 0.0291038, 0.0186934,
783  0.00720046, -0.00537502, -0.0190331, -0.0337736,
784  -0.0495968 },
785  { -0.0340909, -0.0200147, -0.006937, 0.00514208, 0.0162226,
786  0.0263045, 0.0353878, 0.0434725, 0.0505587, 0.0566463,
787  0.0617353, 0.0658257, 0.0689175, 0.0710107, 0.0721054,
788  0.0722014, 0.0712989, 0.0693978, 0.0664981, 0.0625999,
789  0.057703, 0.0518076, 0.0449135, 0.0370209, 0.0281297,
790  0.01824, 0.0073516, -0.00453534, -0.0174209, -0.031305,
791  -0.0461877 },
792  { -0.0369318, -0.0233688, -0.0107221, 0.00100806, 0.0118218,
793  0.0217192, 0.0307001, 0.0387647, 0.0459128, 0.0521444,
794  0.0574597, 0.0618585, 0.0653409, 0.0679069, 0.0695565,
795  0.0702896, 0.0701063, 0.0690066, 0.0669905, 0.0640579,
796  0.0602089, 0.0554435, 0.0497617, 0.0431635, 0.0356488,
797  0.0272177, 0.0178702, 0.0076063, -0.00357405, -0.0156708,
798  -0.028684, -0.0426136 },
799  { -0.038961, -0.025974, -0.0138249, -0.00251362, 0.00795978,
800  0.0175953, 0.026393, 0.0343527, 0.0414747, 0.0477587,
801  0.0532049, 0.0578132, 0.0615836, 0.0645161, 0.0666108,
802  0.0678676, 0.0682866, 0.0678676, 0.0666108, 0.0645161,
803  0.0615836, 0.0578132, 0.0532049, 0.0477587, 0.0414747,
804  0.0343527, 0.026393, 0.0175953, 0.00795978, -0.00251362,
805  -0.0138249, -0.025974, -0.038961 },
806  { -0.0426136, -0.028684, -0.0156708, -0.00357405, 0.0076063,
807  0.0178702, 0.0272177, 0.0356488, 0.0431635, 0.0497617,
808  0.0554435, 0.0602089, 0.0640579, 0.0669905, 0.0690066,
809  0.0701063, 0.0702896, 0.0695565, 0.0679069, 0.0653409,
810  0.0618585, 0.0574597, 0.0521444, 0.0459128, 0.0387647,
811  0.0307001, 0.0217192, 0.0118218, 0.00100806, -0.0107221,
812  -0.0233688, -0.0369318 },
813  { -0.0461877, -0.031305, -0.0174209, -0.00453534, 0.0073516,
814  0.01824, 0.0281297, 0.0370209, 0.0449135, 0.0518076,
815  0.057703, 0.0625999, 0.0664981, 0.0693978, 0.0712989,
816  0.0722014, 0.0721054, 0.0710107, 0.0689175, 0.0658257,
817  0.0617353, 0.0566463, 0.0505587, 0.0434725, 0.0353878,
818  0.0263045, 0.0162226, 0.00514208, -0.006937, -0.0200147,
819  -0.0340909 },
820  { -0.0495968, -0.0337736, -0.0190331, -0.00537502, 0.00720046,
821  0.0186934, 0.0291038, 0.0384316, 0.0466769, 0.0538396,
822  0.0599198, 0.0649174, 0.0688324, 0.0716649, 0.0734149,
823  0.0740823, 0.0736672, 0.0721695, 0.0695892, 0.0659264,
824  0.0611811, 0.0553532, 0.0484427, 0.0404497, 0.0313741,
825  0.021216, 0.00997537, -0.00234785, -0.0157536, -0.0302419 },
826  { -0.0527253, -0.0360083, -0.0204566, -0.00607024, 0.0071508,
827  0.0192065, 0.0300969, 0.039822, 0.0483818, 0.0557763,
828  0.0620054, 0.0670692, 0.0709677, 0.0737009, 0.0752688,
829  0.0756714, 0.0749086, 0.0729806, 0.0698872, 0.0656285,
830  0.0602045, 0.0536151, 0.0458605, 0.0369405, 0.0268552,
831  0.0156046, 0.00318873, -0.0103925, -0.025139 },
832  { -0.0554187, -0.0379037, -0.0216307, -0.00659972, 0.00718917,
833  0.019736, 0.0310408, 0.0411035, 0.0499242, 0.0575028,
834  0.0638394, 0.0689339, 0.0727864, 0.0753968, 0.0767652,
835  0.0768915, 0.0757758, 0.073418, 0.0698181, 0.0649762,
836  0.0588923, 0.0515662, 0.0429982, 0.0331881, 0.0221359,
837  0.00984169, -0.00369458, -0.0184729 },
838  { -0.0574713, -0.0393247, -0.0224833, -0.00694708, 0.0072839,
839  0.0202097, 0.0318302, 0.0421456, 0.0511557, 0.0588607,
840  0.0652604, 0.0703549, 0.0741442, 0.0766284, 0.0778073,
841  0.0776809, 0.0762494, 0.0735127, 0.0694708, 0.0641236,
842  0.0574713, 0.0495137, 0.0402509, 0.029683, 0.0178098,
843  0.00463138, -0.00985222 },
844  { -0.0586081, -0.0400977, -0.0229304, -0.00710623, 0.00737485,
845  0.0205128, 0.0323077, 0.0427595, 0.0518681, 0.0596337,
846  0.0660562, 0.0711355, 0.0748718, 0.077265, 0.078315,
847  0.078022, 0.0763858, 0.0734066, 0.0690842, 0.0634188,
848  0.0564103, 0.0480586, 0.0383639, 0.027326, 0.0149451,
849  0.001221 },
850  { -0.0584615, -0.04, -0.0228763, -0.0070903, 0.00735786,
851  0.0204682, 0.0322408, 0.0426756, 0.0517726, 0.0595318,
852  0.0659532, 0.0710368, 0.0747826, 0.0771906, 0.0782609,
853  0.0779933, 0.076388, 0.0734448, 0.0691639, 0.0635452,
854  0.0565886, 0.0482943, 0.0386622, 0.0276923, 0.0153846 },
855  { -0.0565385, -0.0387458, -0.022215, -0.00694588, 0.00706142,
856  0.0198069, 0.0312907, 0.0415126, 0.0504728, 0.0581712,
857  0.0646078, 0.0697826, 0.0736957, 0.0763469, 0.0777364,
858  0.0778641, 0.07673, 0.0743341, 0.0706765, 0.0657571,
859  0.0595759, 0.0521329, 0.0434281, 0.0334615 },
860  { -0.0521739, -0.0359684, -0.0208357, -0.00677583, 0.00621118,
861  0.0181254, 0.0289667, 0.0387352, 0.0474308, 0.0550536,
862  0.0616036, 0.0670807, 0.071485, 0.0748165, 0.0770751,
863  0.0782609, 0.0783738, 0.0774139, 0.0753811, 0.0722756,
864  0.0680971, 0.0628458, 0.0565217 },
865  { -0.0444664, -0.0311971, -0.0186335, -0.00677583, 0.00437606,
866  0.0148221, 0.0245624, 0.0335968, 0.0419255, 0.0495483,
867  0.0564653, 0.0626765, 0.0681818, 0.0729814, 0.0770751,
868  0.080463, 0.0831451, 0.0851214, 0.0863919, 0.0869565,
869  0.0868154, 0.0859684 },
870  { -0.0321852, -0.0238283, -0.0155309, -0.00729294, 0.000885613,
871  0.00900473, 0.0170644, 0.0250646, 0.0330054, 0.0408868,
872  0.0487087, 0.0564712, 0.0641743, 0.0718179, 0.0794021,
873  0.0869268, 0.0943921, 0.101798, 0.109144, 0.116431,
874  0.123659 },
875  { -0.0136364, -0.0130895, -0.0115174, -0.00892003, -0.00529733,
876  -0.000649351, 0.00502392, 0.0117225, 0.0194463, 0.0281955,
877  0.0379699, 0.0487697, 0.0605947, 0.073445, 0.0873206,
878  0.102221, 0.118148, 0.135099, 0.153076, 0.172078 },
879  { 0.0135338, 0.00200501, -0.00672269, -0.0126493, -0.0157747,
880  -0.0160991, -0.0136223, -0.00834439, -0.000265369, 0.0106148,
881  0.024296, 0.0407784, 0.0600619, 0.0821465, 0.107032,
882  0.134719, 0.165207, 0.198496, 0.234586 },
883  { 0.0526316, 0.0227038, -0.00154799, -0.0201238, -0.0330237,
884  -0.0402477, -0.0417957, -0.0376677, -0.0278638, -0.0123839,
885  0.00877193, 0.0356037, 0.0681115, 0.106295, 0.150155,
886  0.19969, 0.254902, 0.315789 },
887  { 0.108359, 0.0505676, 0.00309598, -0.0340557, -0.0608875,
888  -0.0773994, -0.0835913, -0.0794634, -0.0650155, -0.0402477,
889  -0.00515996, 0.0402477, 0.0959752, 0.162023, 0.23839,
890  0.325077, 0.422085 }
891  };
892 
893  float filtered_baseline[128];
894  float filtered_baseline_derivative[127];
895 
896  // Zero filtered_baseline
897  memset(filtered_baseline, 0, 128 * sizeof(float));
898  // Filter the left edge using (nL, nR) = (0, 16) .. (15, 16) S-G
899  // filters
900  for (size_t i = 0; i < savitzky_golay_n_l_r; i++) {
901  for (size_t j = 0; j < savitzky_golay_n_l_r + 1 + i; j++) {
902  filtered_baseline[i] +=
903  savitzky_golay_coefficient[i][j] * baseline[j];
904  }
905  }
906  // Filter the middle section using the (nL, nR) = (16, 16) S-G
907  // filter, while taking advantage of the symmetry to save 16
908  // multiplications.
909  for (size_t i = savitzky_golay_n_l_r;
910  i < 128 - savitzky_golay_n_l_r; i++) {
911  filtered_baseline[i] =
912  savitzky_golay_coefficient
913  [savitzky_golay_n_l_r][savitzky_golay_n_l_r] * baseline[i];
914  for (size_t j = 0; j < savitzky_golay_n_l_r; j++) {
915  filtered_baseline[i] +=
916  savitzky_golay_coefficient[savitzky_golay_n_l_r][j] *
917  (baseline[i + j - savitzky_golay_n_l_r] +
918  baseline[i - j + savitzky_golay_n_l_r]);
919  }
920 #if 0
921  // Test that the indexing above is correct
922  float test = 0;
923  for (size_t j = 0; j < 2 * savitzky_golay_n_l_r + 1; j++) {
924  test +=
925  savitzky_golay_coefficient[savitzky_golay_n_l_r][j] *
926  baseline[i + j - savitzky_golay_n_l_r];
927  }
928  // test == filtered_baseline[i] should hold now
929 #endif
930  }
931  // Filter the right edge using (nL, nR) = (16, 15) .. (16, 0) S-G
932  // filters
933  for (size_t i = 128 - savitzky_golay_n_l_r; i < 128; i++) {
934  for (size_t j = 0; j < 128 - i + savitzky_golay_n_l_r; j++) {
935  filtered_baseline[i] +=
936  savitzky_golay_coefficient
937  [2 * savitzky_golay_n_l_r + i + 1 - 128][j] *
938  baseline[i + j - savitzky_golay_n_l_r];
939  }
940  }
941  // In lieu of a spearate S-G derivative filter, the finite
942  // difference is used here (since the output is sufficiently
943  // smooth).
944  for (size_t i = 0; i < 127; i++) {
945  filtered_baseline_derivative[i] =
946  filtered_baseline[i + 1] - filtered_baseline[i];
947  }
948 
949  // Calculate the maximum deviation between filtered and unfiltered
950  // baseline, plus the sum square of the derivative.
951 
952  double filtered_baseline_max = 0;
953  double filtered_baseline_derivative_sum_square = 0;
954 
955  for (size_t i = 0; i < 128; i++) {
956  const double d = filtered_baseline[i] - baseline[i];
957 
958  filtered_baseline_max =
959  std::max(filtered_baseline_max,
960  static_cast<double>(fabs(d)));
961  }
962  for (size_t i = 0; i < 127; i++) {
963  filtered_baseline_derivative_sum_square +=
964  filtered_baseline_derivative[i] *
965  filtered_baseline_derivative[i];
966  }
967 
968 #if 0
969  std::cerr << __FILE__ << ':' << __LINE__ << ": "
970  << filtered_baseline_max << ' '
971  << filtered_baseline_derivative_sum_square << std::endl;
972 #endif
973 
974  // Apply the cut
975  return !(filtered_baseline_max >= filteredBaselineMax_ ||
976  filtered_baseline_derivative_sum_square >= filteredBaselineDerivativeSumSquare_);
977 }
978 
979 
980 
981 
982 
983 //Other methods implementation ==============================================
984 //==========================================================================
985 
986 void SiStripAPVRestorer::LoadMeanCMMap(const edm::Event& iEvent){
987  if(useRealMeanCM_){
988  edm::Handle< edm::DetSetVector<SiStripRawDigi> > input;
989  iEvent.getByLabel("siStripDigis","VirginRaw", input);
990  this->CreateCMMapRealPed(*input);
991  } else {
992  edm::Handle< edm::DetSetVector<SiStripProcessedRawDigi> > inputCM;
993  iEvent.getByLabel("MEANAPVCM",inputCM);
994  this->CreateCMMapCMstored(*inputCM);
995  }
996 }
997 
998 
999 void SiStripAPVRestorer::CreateCMMapRealPed(const edm::DetSetVector<SiStripRawDigi>& input){
1000 
1001  MeanCMmap_.erase(MeanCMmap_.begin(), MeanCMmap_.end());
1002 
1003  //std::cout<< "===============================================" << std::endl;
1004 
1005  for ( edm::DetSetVector<SiStripRawDigi>::const_iterator
1006  rawDigis = input.begin(); rawDigis != input.end(); rawDigis++) {
1007  SiStripPedestals::Range detPedestalRange = pedestalHandle->getRange(rawDigis->id);
1008  std::vector<float> MeanCMDetSet;
1009  MeanCMDetSet.clear();
1010 
1011  for(uint16_t APV = 0; APV < rawDigis->size()/128; ++APV){
1012  uint16_t MinPed =0;
1013  for(uint16_t strip = APV*128; strip< (APV+1)*128; ++strip){
1014  uint16_t ped = (uint16_t)pedestalHandle->getPed(strip,detPedestalRange);
1015 
1016  if(ped < MinPed) MinPed = ped;
1017  }
1018  if(MinPed>128) MinPed=128;
1019  MeanCMDetSet.push_back(MinPed);
1020 
1021  }
1022  MeanCMmap_.insert(std::pair<uint32_t, std::vector<float> >(rawDigis->id,MeanCMDetSet));
1023 
1024  }
1025 }
1026 
1027 void SiStripAPVRestorer::CreateCMMapCMstored(const edm::DetSetVector<SiStripProcessedRawDigi>& Input){
1028 
1029  MeanCMmap_.erase(MeanCMmap_.begin(), MeanCMmap_.end());
1030  uint32_t detId;
1031  edm::DetSetVector<SiStripProcessedRawDigi>::const_iterator itInput;
1032  edm::DetSet<SiStripProcessedRawDigi>::const_iterator itCM;
1033  std::vector<float> MeanCMNValue;
1034 
1035  for(itInput = Input.begin(); itInput != Input.end(); ++itInput){
1036  detId = itInput->id;
1037  MeanCMNValue.clear();
1038  for(itCM = itInput->begin(); itCM != itInput->end(); ++itCM) MeanCMNValue.push_back(itCM->adc());
1039  MeanCMmap_.insert(std::pair<uint32_t, std::vector<float> >(detId,MeanCMNValue));
1040  }
1041 }
1042 
1043 std::vector<bool>& SiStripAPVRestorer::GetAPVFlags(){
1044  apvFlagsBool_.clear();
1045  for(size_t i =0; i < apvFlags_.size(); ++i){
1046  if(apvFlags_[i] != "" && !apvFlagsBoolOverride_[i]) apvFlagsBool_.push_back(true);
1047  else apvFlagsBool_.push_back(false);
1048  }
1049  return apvFlagsBool_;
1050 }
1051 
1052 
1053 
1054 
ecalMGPA::adc
int adc(sample_type sample)
get the ADC sample (12 bits)
Definition: EcalMGPASample.h:11
SiStripAPVRestorer::AbnormalBaselineInspect
int16_t AbnormalBaselineInspect(const uint16_t &, std::vector< T > &)
Definition: SiStripAPVRestorer.cc:226
SiStripQuality::IsApvBad
bool IsApvBad(const uint32_t &detid, const short &apvNb) const
Definition: SiStripQuality.cc:617
SiStripAPVRestorer::BaselineFollowerInspect
int16_t BaselineFollowerInspect(const uint16_t &, std::vector< T > &)
Definition: SiStripAPVRestorer.cc:142
SiStripAPVRestorer::CreateCMMapRealPed
void CreateCMMapRealPed(const edm::DetSetVector< SiStripRawDigi > &)
Definition: SiStripAPVRestorer.cc:999
SiStripAPVRestorer::restore
void restore(const uint16_t &, std::vector< int16_t > &)
Definition: SiStripAPVRestorer.cc:111
SiStripAPVRestorer::init
void init(const edm::EventSetup &es)
Definition: SiStripAPVRestorer.cc:45
mps_fire.i
i
Definition: mps_fire.py:156
SiStripAPVRestorer::GetAPVFlags
std::vector< bool > & GetAPVFlags()
Definition: SiStripAPVRestorer.cc:1043
SiStripAPVRestorer::FlatRegionsFinder
bool FlatRegionsFinder(const std::vector< int16_t > &, DigiMap &, const uint16_t &)
Definition: SiStripAPVRestorer.cc:375
SiStripAPVRestorer::LoadMeanCMMap
void LoadMeanCMMap(const edm::Event &)
Definition: SiStripAPVRestorer.cc:986
hpstanc_transforms.max
max
Definition: hpstanc_transforms.py:5
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
SiStripQuality::IsStripBad
bool IsStripBad(const uint32_t &detid, const short &strip) const
Definition: SiStripQuality.cc:625
Input
#define Input(cl)
Definition: vmac.h:188
MessageLogger_cfi.cerr
cerr
Definition: MessageLogger_cfi.py:547
DigiMap
std::map< uint16_t, int16_t > DigiMap
Definition: SiStripAPVRestorer.h:25
SiStripAPVRestorer::filteredBaselineDerivativeSumSquare_
double filteredBaselineDerivativeSumSquare_
Definition: SiStripAPVRestorer.h:138
SiStripPedestals::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripPedestals.h:52
SiStripAPVRestorer::pedestalHandle
edm::ESHandle< SiStripPedestals > pedestalHandle
Definition: SiStripAPVRestorer.h:90
SiStripAPVRestorer::Cleaner_MonotonyChecker
void Cleaner_MonotonyChecker(DigiMap &)
Definition: SiStripAPVRestorer.cc:496
SiStripAPVRestorer::SmoothedMaps_
std::map< uint16_t, DigiMap > SmoothedMaps_
Definition: SiStripAPVRestorer.h:98
input
static std::string const input
Definition: EdmProvDump.cc:44
DigiMapIter
std::map< uint16_t, int16_t >::iterator DigiMapIter
Definition: SiStripAPVRestorer.h:27
SiStripPedestals::getPed
float getPed(const uint16_t &strip, const Range &range) const
Definition: SiStripPedestals.cc:51
SiStripAPVRestorer::apvFlags_
std::vector< std::string > apvFlags_
Definition: SiStripAPVRestorer.h:93
SiStripNoises::getNoiseFast
static float getNoiseFast(const uint16_t &strip, const Range &range)
Definition: SiStripNoises.h:67
SiStripAPVRestorer::BaselineFollowerRestore
void BaselineFollowerRestore(const uint16_t &, const uint16_t &, const float &, std::vector< int16_t > &)
Definition: SiStripAPVRestorer.cc:308
SiStripAPVRestorer::CheckBaseline
bool CheckBaseline(const std::vector< int16_t > &) const
Definition: SiStripAPVRestorer.cc:688
iEvent
int iEvent
Definition: GenABIO.cc:230
SiStripAPVRestorer::qualityHandle
edm::ESHandle< SiStripQuality > qualityHandle
Definition: SiStripAPVRestorer.h:84
SiStripAPVRestorer::BaselineFollower
void BaselineFollower(DigiMap &, std::vector< int16_t > &, const float &)
Definition: SiStripAPVRestorer.cc:632
SiStripAPVRestorer::InspectAndRestore
int16_t InspectAndRestore(const uint32_t &, const uint16_t &, std::vector< int16_t > &, std::vector< int16_t > &, const std::vector< std::pair< short, float > > &)
Definition: SiStripAPVRestorer.cc:73
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SiStripAPVRestorer::NullInspect
int16_t NullInspect(const uint16_t &, std::vector< T > &)
Definition: SiStripAPVRestorer.cc:270
min
T min(T a, T b)
Definition: MathUtil.h:58
digi_MixPreMix_cfi.strip
strip
Definition: digi_MixPreMix_cfi.py:15
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:416
edm::DetSetVector::end
iterator end()
Return the off-the-end iterator.
Definition: DetSetVector.h:361
SiStripAPVRestorer::BaselineAndSaturationInspect
int16_t BaselineAndSaturationInspect(const uint16_t &, std::vector< T > &)
Definition: SiStripAPVRestorer.cc:186
SiStripAPVRestorer::FlatRestore
void FlatRestore(const uint16_t &, const uint16_t &, std::vector< int16_t > &)
Definition: SiStripAPVRestorer.cc:354
SiStripAPVRestorer::Cleaner_HighSlopeChecker
void Cleaner_HighSlopeChecker(DigiMap &)
Definition: SiStripAPVRestorer.cc:600
SiStripAPVRestorer::apvFlagsBoolOverride_
std::vector< bool > apvFlagsBoolOverride_
Definition: SiStripAPVRestorer.h:95
SiStripAPVRestorer::noiseHandle
edm::ESHandle< SiStripNoises > noiseHandle
Definition: SiStripAPVRestorer.h:87
SiStripAPVRestorer::inspect
int16_t inspect(const uint32_t &, const uint16_t &, std::vector< int16_t > &, const std::vector< std::pair< short, float > > &)
Definition: SiStripAPVRestorer.cc:80
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:56
SiStripAPVRestorer::BaselineCleaner
void BaselineCleaner(const std::vector< int16_t > &, DigiMap &, const uint16_t &)
Definition: SiStripAPVRestorer.cc:474
SiStripAPVRestorer::badAPVs_
std::vector< bool > badAPVs_
Definition: SiStripAPVRestorer.h:97
SiStripAPVRestorer::apvFlagsBool_
std::vector< bool > apvFlagsBool_
Definition: SiStripAPVRestorer.h:94
SiStripBadStrip::getRange
const Range getRange(const uint32_t detID) const
Definition: SiStripBadStrip.cc:26
SiStripNoises::getRange
const Range getRange(const uint32_t detID) const
Definition: SiStripNoises.cc:34
SiStripBadStrip::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripBadStrip.h:54
SiStripAPVRestorer::median_
std::vector< float > median_
Definition: SiStripAPVRestorer.h:96
SiStripAPVRestorer::CreateCMMapCMstored
void CreateCMMapCMstored(const edm::DetSetVector< SiStripProcessedRawDigi > &)
Definition: SiStripAPVRestorer.cc:1027
edm::DetSetVector::begin
iterator begin()
Return an iterator to the first DetSet.
Definition: DetSetVector.h:346
SiStripPedestals::getRange
const Range getRange(const uint32_t &detID) const
Definition: SiStripPedestals.cc:27
SiStripAPVRestorer::Cleaner_LocalMinimumAdder
void Cleaner_LocalMinimumAdder(const std::vector< int16_t > &, DigiMap &, const uint16_t &)
Definition: SiStripAPVRestorer.cc:522
SiStripNoises::Range
std::pair< ContainerIterator, ContainerIterator > Range
Definition: SiStripNoises.h:48
edm::DetSet::const_iterator
collection_type::const_iterator const_iterator
Definition: DetSet.h:33
edm::DetSetVector::const_iterator
collection_type::const_iterator const_iterator
Definition: DetSetVector.h:104
SiStripAPVRestorer::SiStripAPVRestorer
SiStripAPVRestorer(const edm::ParameterSet &conf)
Definition: SiStripAPVRestorer.cc:9
Generated for CMSSW Reference Manual by doxygen 1.8.11