CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoLocalTracker/SiStripZeroSuppression/src/SiStripAPVRestorer.cc

Go to the documentation of this file.
00001 #include "RecoLocalTracker/SiStripZeroSuppression/interface/SiStripAPVRestorer.h"
00002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00003 
00004 #include <cmath>
00005 #include <iostream>
00006 #include <algorithm>
00007 
00008 
00009 SiStripAPVRestorer::SiStripAPVRestorer(const edm::ParameterSet& conf):
00010   quality_cache_id(-1), noise_cache_id(-1), pedestal_cache_id(-1),
00011   ForceNoRestore_(conf.getParameter<bool>("ForceNoRestore")),
00012   SelfSelectRestoreAlgo_(conf.getParameter<bool>("SelfSelectRestoreAlgo")),
00013   InspectAlgo_(conf.getParameter<std::string>("APVInspectMode")),
00014   RestoreAlgo_(conf.getParameter<std::string>("APVRestoreMode")),
00015   useRealMeanCM_(conf.getParameter<bool>("useRealMeanCM")),
00016   fraction_(conf.getParameter<double>("Fraction")),
00017   deviation_(conf.getParameter<uint32_t>("Deviation")),
00018   restoreThreshold_(conf.getParameter<double>("restoreThreshold")),
00019   DeltaCMThreshold_(conf.getParameter<uint32_t>("DeltaCMThreshold")),
00020   nSigmaNoiseDerTh_(conf.getParameter<uint32_t>("nSigmaNoiseDerTh")),
00021   consecThreshold_(conf.getParameter<uint32_t>("consecThreshold")),
00022   hitStripThreshold_(conf.getParameter<uint32_t>("hitStripThreshold")),
00023   nSmooth_(conf.getParameter<uint32_t>("nSmooth")),
00024   minStripsToFit_(conf.getParameter<uint32_t>("minStripsToFit")),
00025   distortionThreshold_(conf.getParameter<uint32_t>("distortionThreshold")),
00026   CutToAvoidSignal_(conf.getParameter<double>("CutToAvoidSignal")),
00027   nSaturatedStrip_(conf.getParameter<uint32_t>("nSaturatedStrip")),
00028   ApplyBaselineCleaner_(conf.getParameter<bool>("ApplyBaselineCleaner")),
00029   slopeX_(conf.getParameter<int32_t>("slopeX")),
00030   slopeY_(conf.getParameter<int32_t>("slopeY")),
00031   CleaningSequence_(conf.getParameter<uint32_t>("CleaningSequence")),
00032   ApplyBaselineRejection_(conf.getParameter<bool>("ApplyBaselineRejection")),
00033   MeanCM_(conf.getParameter<int32_t>("MeanCM")),
00034   filteredBaselineMax_(conf.getParameter<double>("filteredBaselineMax")),
00035   filteredBaselineDerivativeSumSquare_(conf.getParameter<double>("filteredBaselineDerivativeSumSquare"))  
00036 
00037 {
00038   apvFlags_.clear();
00039   median_.clear();
00040   SmoothedMaps_.clear();
00041   BaselineMap_.erase(BaselineMap_.begin(), BaselineMap_.end());
00042 }
00043 
00044 
00045 void SiStripAPVRestorer::init(const edm::EventSetup& es){
00046   uint32_t n_cache_id = es.get<SiStripNoisesRcd>().cacheIdentifier();
00047   uint32_t q_cache_id = es.get<SiStripQualityRcd>().cacheIdentifier();
00048   uint32_t p_cache_id = es.get<SiStripPedestalsRcd>().cacheIdentifier();
00049   
00050   if(n_cache_id != noise_cache_id) {
00051     es.get<SiStripNoisesRcd>().get( noiseHandle );
00052     noise_cache_id = n_cache_id;
00053   } else {
00054     noise_cache_id = n_cache_id;
00055   }
00056   if(q_cache_id != quality_cache_id) {
00057     es.get<SiStripQualityRcd>().get( qualityHandle );
00058     quality_cache_id = q_cache_id;
00059   }else {
00060     quality_cache_id = q_cache_id;
00061   }
00062   
00063   if(p_cache_id != pedestal_cache_id) {
00064                 es.get<SiStripPedestalsRcd>().get( pedestalHandle );
00065                 pedestal_cache_id = p_cache_id;
00066   }else {
00067     pedestal_cache_id = p_cache_id;
00068   }
00069   
00070 }
00071 
00072  
00073 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 ){
00074   int16_t nAPVFlagged = this->inspect(detId, firstAPV, rawDigisPedSubtracted, vmedians);
00075   this->restore(firstAPV, processedRawDigi);
00076   return nAPVFlagged;
00077 }
00078 
00079 
00080 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) {
00081   
00082   detId_ = detId;
00083   
00084   apvFlagsBool_.clear();
00085   apvFlagsBoolOverride_.clear();
00086   apvFlagsBoolOverride_.insert(apvFlagsBoolOverride_.begin(), 6, false);
00087   apvFlags_.clear();
00088   apvFlags_.insert(apvFlags_.begin(), 6, "");
00089   median_.clear();
00090   median_.insert(median_.begin(), 6, -999);
00091   badAPVs_.clear();
00092   badAPVs_.insert(badAPVs_.begin(), 6, false);
00093   SmoothedMaps_.erase(SmoothedMaps_.begin(),  SmoothedMaps_.end());
00094   BaselineMap_.erase(BaselineMap_.begin(), BaselineMap_.end()); 
00095     
00096   for(size_t i=0; i< vmedians.size(); ++i){
00097          short APV =  vmedians[i].first;
00098          median_[APV]= vmedians[i].second;
00099          badAPVs_[APV] = qualityHandle->IsApvBad(detId_, APV);
00100   }
00101         
00102   if(InspectAlgo_=="BaselineFollower") return this->BaselineFollowerInspect(firstAPV, digis); 
00103   if(InspectAlgo_=="AbnormalBaseline") return this->AbnormalBaselineInspect(firstAPV, digis);
00104   if(InspectAlgo_=="Null") return this->NullInspect(firstAPV, digis);
00105   if(InspectAlgo_=="BaselineAndSaturation") return this->BaselineAndSaturationInspect(firstAPV, digis);
00106   throw cms::Exception("Unregistered Inspect Algorithm") << "SiStripAPVRestorer possibilities: (Null), (AbnormalBaseline),(BaselineFollower)";
00107   
00108 }
00109 
00110 
00111 void SiStripAPVRestorer::restore(const uint16_t& firstAPV, std::vector<int16_t>& digis ) {
00112         
00113   if(ForceNoRestore_) return;
00114   
00115   for( uint16_t APV=firstAPV; APV< digis.size()/128 + firstAPV; ++APV){
00116     std::string algoToUse = *( apvFlags_.begin() + APV );
00117     
00118     if ( algoToUse != ""){
00119       if(!SelfSelectRestoreAlgo_) algoToUse = RestoreAlgo_;
00120  
00121       if(algoToUse=="Flat"){
00122         this->FlatRestore(APV, firstAPV, digis);
00123       }else if(algoToUse=="BaselineFollower"){
00124         this->BaselineFollowerRestore(APV, firstAPV, median_[APV], digis);
00125       }else{
00126         throw cms::Exception("Unregistered Restore Algorithm") << "SiStripAPVRestorer possibilities: (Flat), (BaselineFollower)";
00127       }
00128       
00129       
00130     }
00131   }
00132   
00133 }
00134 
00135 
00136 //Inspect method implementation ========================================================================================================================================================================
00137 //======================================================================================================================================================================================================
00138 //======================================================================================================================================================================================================
00139 
00140 template<typename T>
00141 inline
00142 int16_t SiStripAPVRestorer::BaselineFollowerInspect(const uint16_t& firstAPV, std::vector<T>& digis){
00143   std::vector<T> singleAPVdigi;  
00144   int16_t nAPVflagged = 0;
00145   
00146   CMMap::iterator itCMMap;
00147   if(useRealMeanCM_) itCMMap = MeanCMmap_.find(detId_);
00148   
00149   for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
00150 
00151     DigiMap smoothedmap;
00152     smoothedmap.erase(smoothedmap.begin(), smoothedmap.end());
00153 
00154     if(!badAPVs_[APV]){
00155       float MeanAPVCM = MeanCM_;
00156       if(useRealMeanCM_&&itCMMap!= MeanCMmap_.end()) MeanAPVCM =(itCMMap->second)[APV];
00157     
00158       singleAPVdigi.clear(); 
00159       for(int16_t strip = (APV-firstAPV)*128; strip < (APV-firstAPV+1)*128; ++strip){
00160         singleAPVdigi.push_back(digis[strip]); 
00161       }
00162    
00163    
00164       float DeltaCM = median_[APV] - MeanAPVCM; 
00165       
00166       //std::cout << "Delta CM: " << DeltaCM << " CM: " << median_[APV] << " detId " << (uint32_t) detId_ << std::endl;         
00167       if(DeltaCM < 0 && std::abs(DeltaCM) > DeltaCMThreshold_){
00168       
00169         bool isFlat = FlatRegionsFinder(singleAPVdigi,smoothedmap,APV);
00170         if(!isFlat){
00171               apvFlags_[APV]= "BaselineFollower";    //specify any algo to make the restore
00172               nAPVflagged++;
00173         }
00174       } 
00175       
00176     } 
00177     SmoothedMaps_.insert(SmoothedMaps_.end(), std::pair<uint16_t, DigiMap>(APV, smoothedmap));
00178    }
00179   
00180   return nAPVflagged;
00181 }
00182 
00183 //======================================================================================================================================================================================================
00184 template<typename T>
00185 inline
00186 int16_t SiStripAPVRestorer::BaselineAndSaturationInspect(const uint16_t& firstAPV, std::vector<T>& digis){
00187   std::vector<T> singleAPVdigi;
00188   singleAPVdigi.clear();
00189   
00190   
00191   int16_t nAPVflagged = 0;
00192   
00193   CMMap::iterator itCMMap;
00194   if(useRealMeanCM_) itCMMap = MeanCMmap_.find(detId_);
00195   
00196   for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
00197      apvFlags_.push_back( "" );
00198     if(!badAPVs_[APV]){
00199      float MeanAPVCM = MeanCM_;
00200      if(useRealMeanCM_&&itCMMap!= MeanCMmap_.end()) MeanAPVCM =(itCMMap->second)[APV];
00201     
00202      singleAPVdigi.clear();
00203    
00204      uint16_t nSatStrip =0;
00205      for(int16_t strip = (APV-firstAPV)*128; strip < (APV-firstAPV+1)*128; ++strip){
00206         singleAPVdigi.push_back(digis[strip]); 
00207         if(digis[strip] >=1023) ++nSatStrip;
00208       }
00209          
00210      float DeltaCM = median_[APV] -MeanAPVCM; 
00211     
00212     
00213      if(DeltaCM < 0 && std::abs(DeltaCM) > DeltaCMThreshold_&&nSatStrip>= nSaturatedStrip_){
00214        apvFlags_[APV] = RestoreAlgo_;    //specify any algo to make the restore
00215        nAPVflagged++;
00216      } 
00217     }   
00218   }
00219   
00220   return nAPVflagged;
00221 }
00222 
00223 //======================================================================================================================================================================================================
00224 template<typename T>
00225 inline
00226 int16_t SiStripAPVRestorer::AbnormalBaselineInspect( const uint16_t& firstAPV, std::vector<T>& digis){
00227 
00228   SiStripQuality::Range detQualityRange = qualityHandle->getRange(detId_);
00229   
00230   typename std::vector<T>::iterator fs;
00231   
00232   int16_t nAPVflagged=0;
00233   
00234   CMMap::iterator itCMMap;
00235   if(useRealMeanCM_) itCMMap = MeanCMmap_.find(detId_);
00236   
00237   
00238   int devCount = 0, qualityCount = 0, minstrip = 0; 
00239  for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
00240     apvFlags_.push_back( "" );
00241     if(!badAPVs_[APV]){
00242       float MeanAPVCM = MeanCM_;
00243       if(useRealMeanCM_&&itCMMap!= MeanCMmap_.end()) MeanAPVCM =(itCMMap->second)[APV];
00244       for (uint16_t istrip=APV*128; istrip<(APV+1)*128; ++istrip){
00245         fs = digis.begin() + istrip-firstAPV*128;
00246         if ( !qualityHandle->IsStripBad(detQualityRange,istrip) ){
00247                qualityCount++; 
00248                if ( std::abs((int) *fs - MeanAPVCM) > (int)deviation_ ){ 
00249                 devCount++;
00250                     minstrip = std::min((int) *fs, minstrip);
00251            }
00252          }
00253       }
00254     
00255       if( devCount > fraction_ * qualityCount ) {
00256         apvFlags_[APV] = RestoreAlgo_;      //specify any algo to make the restore
00257         nAPVflagged++;
00258       } 
00259     } 
00260   }
00261   
00262   return nAPVflagged;
00263   
00264 }
00265 
00266 
00267 //======================================================================================================================================================================================================
00268 template<typename T>
00269 inline
00270 int16_t SiStripAPVRestorer::NullInspect(const uint16_t& firstAPV, std::vector<T>& digis){
00271 
00272   SiStripQuality::Range detQualityRange = qualityHandle->getRange(detId_);
00273 
00274   typename std::vector<T>::iterator fs;
00275 
00276   int16_t nAPVflagged = 0;
00277 
00278   for(uint16_t APV=firstAPV ; APV< digis.size()/128 + firstAPV; ++APV){
00279    apvFlags_.push_back( "" );
00280    if(!badAPVs_[APV]){ 
00281      int zeroCount = 0, qualityCount = 0; 
00282      for (uint16_t istrip=APV*128; istrip<(APV+1)*128; ++istrip){
00283        fs = digis.begin() + istrip-firstAPV*128;
00284        if ( !qualityHandle->IsStripBad(detQualityRange,istrip) ){
00285         qualityCount++; 
00286         if ( (int) *fs < 1 ) zeroCount++;
00287        }
00288       }
00289     
00290       if( zeroCount > restoreThreshold_ * qualityCount ) {
00291         apvFlags_[APV] = RestoreAlgo_;     //specify any algo to make the restore
00292         nAPVflagged++;
00293       } 
00294    } 
00295    }
00296  
00297   return nAPVflagged;
00298 
00299 }
00300 
00301 
00302 
00303 //Restore method implementation ========================================================================================================================================================================
00304 //======================================================================================================================================================================================================
00305 //======================================================================================================================================================================================================
00306 
00307 inline
00308 void SiStripAPVRestorer::BaselineFollowerRestore(const uint16_t& APVn, const uint16_t& firstAPV, const float& median, std::vector<int16_t>& digis){
00309   //typename std::vector<T>::iterator firstStrip(digis.begin() + APVn*128), lastStrip(firstStrip + 128), actualStrip;
00310   
00311   
00312   std::vector<int16_t> baseline;
00313   baseline.clear();
00314   baseline.insert(baseline.begin(),128, 0);
00315         
00316          
00317  
00318     
00319   //============================= Find Flat Regions & Interpolating the baseline & subtracting the baseline  =================  
00320   
00321   if(SmoothedMaps_.size()){
00322     std::map<uint16_t, DigiMap >::iterator itSmootedMap = SmoothedMaps_.find(APVn);     
00323     this->BaselineFollower(itSmootedMap->second, baseline, median);
00324   } else {
00325     //median=0;
00326     DigiMap  smoothedpoints;
00327     std::vector<int16_t> singleAPVdigi;
00328     singleAPVdigi.clear(); 
00329     for(int16_t strip = (APVn-firstAPV)*128; strip < (APVn-firstAPV+1)*128; ++strip) singleAPVdigi.push_back(digis[strip]); 
00330     this->FlatRegionsFinder(singleAPVdigi,smoothedpoints, APVn);
00331     this->BaselineFollower(smoothedpoints, baseline, median);           
00332     
00333   }     
00334 
00335   if(ApplyBaselineRejection_){
00336     if(CheckBaseline(baseline)) apvFlagsBoolOverride_[APVn] = true;
00337   }
00338   
00339   //============================= subtracting the baseline =============================================
00340   
00341   for(int16_t itStrip= 0 ; itStrip< 128; ++itStrip){
00342     digis[(APVn-firstAPV)*128+itStrip] -= baseline[itStrip] - median;
00343   }
00344   
00345                 
00346   //============================= storing baseline to the map =============================================     
00347   BaselineMap_.insert(BaselineMap_.end(),  std::pair< uint16_t, std::vector < int16_t> >(APVn, baseline));
00348   
00349 }
00350 
00351 
00352 //======================================================================================================================================================================================================
00353 inline
00354 void SiStripAPVRestorer::FlatRestore(const uint16_t& APVn, const uint16_t& firstAPV, std::vector<int16_t>& digis ){
00355  
00356   std::vector<int16_t> baseline;
00357   baseline.clear();
00358   baseline.insert(baseline.begin(),128, 150);
00359   baseline[0]=0; baseline[127]=0;
00360   BaselineMap_.insert(BaselineMap_.end(),  std::pair< uint16_t, std::vector < int16_t> >(APVn, baseline));  
00361   
00362   for(int16_t itStrip= 0 ; itStrip< 128; ++itStrip){
00363     digis[(APVn-firstAPV)*128+itStrip] = baseline[itStrip];
00364   }
00365  
00366   
00367 }
00368 
00369 
00370 
00371 //Baseline calculation implementation ==================================================================================================================================================================
00372 //======================================================================================================================================================================================================
00373 //======================================================================================================================================================================================================
00374 
00375 bool inline SiStripAPVRestorer::FlatRegionsFinder(const std::vector<int16_t>& adcs, DigiMap& smoothedpoints, const uint16_t& APVn){
00376   SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detId_);
00377   
00378   DigiMap consecpoints;
00379   DigiMapIter itConsecpoints, itSmoothedpoints;
00380   consecpoints.erase(consecpoints.begin(), consecpoints.end());
00381   smoothedpoints.erase(smoothedpoints.begin(), smoothedpoints.end());
00382   
00383    
00384   //============================= Height above local minimum ===============================                    
00385   std::vector<float> adcsLocalMinSubtracted;
00386   adcsLocalMinSubtracted.clear();
00387   adcsLocalMinSubtracted.insert(adcsLocalMinSubtracted.begin(), 128,0);
00388   for(uint32_t istrip=0; istrip<128; ++istrip) {
00389     float localmin = 999.9;             
00390     for(uint16_t jstrip=std::max(0,(int)(istrip-nSmooth_/2)); jstrip<std::min(128,(int)(istrip+nSmooth_/2)); ++jstrip) {
00391       float nextvalue = adcs[jstrip];
00392       if(nextvalue < localmin) localmin=nextvalue;                      
00393     }
00394     adcsLocalMinSubtracted[istrip] = adcs[istrip] - localmin;
00395   }
00396   
00397   
00398   //============================= Find regions with stable slopes ========================
00399   std::vector<uint16_t> nConsStrip;
00400   nConsStrip.clear();
00401   
00402   //Creating maps with all the neighborhood strip and putting in a nCosntStip vector how many we have
00403   uint16_t consecStrips=0;
00404   for(uint32_t istrip=0; istrip<128; ++istrip) {    
00405     int16_t adc = adcs[istrip]; 
00406  
00407    //if( adcsLocalMinSubtracted[istrip] < nSigmaNoiseDerTh_ * (float)noiseHandle->getNoise(istrip+APVn*128,detNoiseRange) && (adc - median) < hitStripThreshold_){
00408    if( adcsLocalMinSubtracted[istrip] < nSigmaNoiseDerTh_ * (float)noiseHandle->getNoiseFast(istrip+APVn*128,detNoiseRange)){
00409       consecpoints.insert(consecpoints.end(), std::pair<uint16_t, int16_t >(istrip, adc));
00410       ++consecStrips;
00411     }else if (consecStrips >0){
00412       nConsStrip.push_back(consecStrips);
00413       consecStrips = 0;
00414     }    
00415   }                     
00416 
00417   //to cope with the last flat region of the APV
00418   if(consecStrips >0) nConsStrip.push_back(consecStrips);
00419 
00420   //removing from the map the fist and last points in wide flat regions and erasing from the map too small regions
00421   itConsecpoints = consecpoints.begin();
00422   float MinSmoothValue=20000., MaxSmoothValue=0.;
00423   for(std::vector<uint16_t>::iterator itnConsStrip = nConsStrip.begin(); itnConsStrip < nConsStrip.end(); ++itnConsStrip){
00424     
00425     consecStrips = *itnConsStrip;
00426     if(consecStrips >=consecThreshold_){
00427       ++itConsecpoints;  //skipping first point
00428       uint16_t nFirstStrip = itConsecpoints->first;
00429       uint16_t nLastStrip;
00430       float smoothValue = 0.0;
00431       float stripCount =1;
00432       for(uint16_t n =0; n < consecStrips-2; ++n){
00433                 smoothValue += itConsecpoints->second;
00434                 if(stripCount == consecThreshold_){
00435                   smoothValue /= (float)stripCount;
00436                   nLastStrip = nFirstStrip + stripCount -1;                                                 
00437                   smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nFirstStrip, smoothValue));
00438                   smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nLastStrip, smoothValue));
00439                   if(smoothValue > MaxSmoothValue) MaxSmoothValue = smoothValue;
00440                   if(smoothValue < MinSmoothValue) MinSmoothValue = smoothValue;
00441                   nFirstStrip = nLastStrip+1;
00442                   smoothValue=0;
00443                   stripCount=0;
00444                 }
00445                 ++stripCount;
00446                 ++itConsecpoints;
00447      }
00448      ++itConsecpoints;  //and putting the pointer to the new seies of point 
00449       
00450      if(stripCount>1) {
00451      //if(smoothValue>0){
00452                 --stripCount;
00453                 smoothValue /= (float)(stripCount);
00454                 nLastStrip = nFirstStrip + stripCount -1;
00455                 smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nFirstStrip, smoothValue));
00456                 smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(nLastStrip, smoothValue));
00457                 if(smoothValue > MaxSmoothValue) MaxSmoothValue = smoothValue;
00458                 if(smoothValue < MinSmoothValue) MinSmoothValue = smoothValue;
00459      }
00460    } else{
00461       for(int n =0; n< consecStrips ; ++n) ++itConsecpoints;
00462    }
00463   }
00464   
00465         
00466   if( (MaxSmoothValue-MinSmoothValue) > distortionThreshold_){
00467         if(ApplyBaselineCleaner_) this->BaselineCleaner(adcs, smoothedpoints, APVn);
00468         return false;
00469   }
00470   return true;
00471 }
00472 
00473 
00474 void inline SiStripAPVRestorer::BaselineCleaner(const std::vector<int16_t>& adcs, DigiMap& smoothedpoints, const uint16_t& APVn){
00475   
00476   if(CleaningSequence_==0) {  //default sequence used up to now
00477     this->Cleaner_HighSlopeChecker(smoothedpoints);
00478     this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
00479   }else if(CleaningSequence_==1){     
00480     this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
00481     this->Cleaner_HighSlopeChecker(smoothedpoints);
00482     this->Cleaner_MonotonyChecker(smoothedpoints);
00483   }else if(CleaningSequence_==2){
00484     this->Cleaner_HighSlopeChecker(smoothedpoints);
00485   }else if(CleaningSequence_==3){
00486     this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
00487     this->Cleaner_HighSlopeChecker(smoothedpoints);    
00488   }else{
00489     this->Cleaner_HighSlopeChecker(smoothedpoints);
00490     this->Cleaner_LocalMinimumAdder(adcs, smoothedpoints, APVn);
00491   }     
00492     
00493 }
00494 
00495 
00496 void inline SiStripAPVRestorer::Cleaner_MonotonyChecker(DigiMap& smoothedpoints){
00497 //Removing points without monotony
00498         //--------------------------------------------------------------------------------------------------
00499          if(smoothedpoints.size() < 3) return;         
00500          DigiMapIter itSmoothedpoints, itSmoothedpointsNext, itSmoothedpointsNextNext, itSmoothedpointsBegin, itSmoothedpointsEnd;
00501          
00502         itSmoothedpoints=smoothedpoints.begin();
00503         while (smoothedpoints.size() > 3 && itSmoothedpoints != --(--(smoothedpoints.end()))) { //while we are not at the last point
00504                 // get info about current and next points
00505                 itSmoothedpointsNext = itSmoothedpoints;
00506                 ++itSmoothedpointsNext;
00507                 itSmoothedpointsNextNext = itSmoothedpointsNext;
00508                 ++itSmoothedpointsNextNext;
00509                 float adc1 = itSmoothedpoints->second;
00510                 float adc2 = itSmoothedpointsNext->second;
00511                 float adc3 = itSmoothedpointsNextNext->second;
00512                 
00513                 if((adc2-adc1) > hitStripThreshold_ && (adc2-adc3) > hitStripThreshold_){
00514                   smoothedpoints.erase(itSmoothedpointsNext);
00515                 }else {
00516                   ++itSmoothedpoints;
00517                 }
00518                         
00519         }
00520 }
00521 
00522 void inline SiStripAPVRestorer::Cleaner_LocalMinimumAdder(const std::vector<int16_t>& adcs, DigiMap& smoothedpoints, const uint16_t& APVn){
00523   SiStripNoises::Range detNoiseRange = noiseHandle->getRange(detId_);
00524   //inserting extra point is case of local minimum
00525         //--------------------------------------------------------------------------------------------------
00526         // these should be reset now for the point-insertion that follows
00527         
00528         DigiMapIter itSmoothedpoints, itSmoothedpointsNext, itSmoothedpointsBegin, itSmoothedpointsEnd; 
00529         if(smoothedpoints.size() >= 2){
00530  
00531         itSmoothedpointsBegin = smoothedpoints.begin();
00532         itSmoothedpointsEnd = --(smoothedpoints.end());
00533                 for(itSmoothedpoints = itSmoothedpointsBegin; itSmoothedpoints != itSmoothedpointsEnd; ++itSmoothedpoints){  
00534                 itSmoothedpointsNext = itSmoothedpoints;
00535                 ++itSmoothedpointsNext;
00536                 float strip1 = itSmoothedpoints->first;
00537                 float strip2 = itSmoothedpointsNext->first;
00538                 float adc1 = itSmoothedpoints->second;
00539                 float adc2 = itSmoothedpointsNext->second;
00540                 float m = (adc2 -adc1)/(strip2 -strip1);
00541     
00542                 //2,4
00543                 if((strip2 - strip1) >slopeX_ && abs(adc1 -adc2) >slopeY_){
00544                         float itStrip = 1;
00545                         float strip = itStrip + strip1;
00546                         while(strip < strip2){
00547                                 
00548                         float adc = adcs[strip];
00549                         if( adc < (adc1 + m * itStrip - 2 * (float)noiseHandle->getNoiseFast(strip+APVn*128,detNoiseRange))){
00550                                                 //std::cout << "applying correction strip: " << strip + APVn*128 << " adc " << adc << " detId: " << detId_ << std::endl;
00551                                                 smoothedpoints.insert(itSmoothedpointsNext, std::pair<uint16_t, int16_t >(strip,adc));
00552                                                 ++itSmoothedpoints;
00553                                                 ++itSmoothedpointsNext;
00554                                                 itSmoothedpointsEnd = --(smoothedpoints.end());
00555                                         } 
00556                                         ++itStrip;
00557                                         ++strip;
00558                                 }
00559                         
00560 
00561                 }
00562                 }
00563         }
00564         
00565         
00566     itSmoothedpointsBegin = smoothedpoints.begin();
00567     itSmoothedpointsEnd = --(smoothedpoints.end());
00568     uint16_t firstStripFlat = itSmoothedpointsBegin->first;
00569     uint16_t lastStripFlat = itSmoothedpointsEnd->first;
00570     int16_t firstStripFlatADC= itSmoothedpointsBegin->second;
00571     int16_t lastStripFlatADC= itSmoothedpointsEnd->second;
00572         
00573     itSmoothedpoints = itSmoothedpointsBegin;
00574     if(firstStripFlat >3){
00575                 float strip = 0;
00576         while(strip < firstStripFlat){
00577                         float adc = adcs[strip];
00578             if( adc < ( firstStripFlatADC - 2 * (float)noiseHandle->getNoiseFast(strip+APVn*128,detNoiseRange))){
00579                                         smoothedpoints.insert(itSmoothedpoints, std::pair<uint16_t, int16_t >(strip,adc));
00580                                         ++itSmoothedpoints;
00581                         } 
00582                         ++strip;
00583                 }
00584         }
00585         
00586         itSmoothedpoints = itSmoothedpointsEnd;
00587         if(lastStripFlat <125){
00588                 float strip = lastStripFlat+1;
00589         while(strip < 128){
00590                         float adc = adcs[strip];
00591             if( adc < ( lastStripFlatADC - 2 * (float)noiseHandle->getNoiseFast(strip+APVn*128,detNoiseRange))){
00592                 smoothedpoints.insert(smoothedpoints.end(), std::pair<uint16_t, int16_t >(strip,adc));
00593                         } 
00594                         ++strip;
00595                 }
00596         }
00597 }
00598 
00599 
00600 void inline SiStripAPVRestorer::Cleaner_HighSlopeChecker(DigiMap& smoothedpoints){
00601   //Removing points in the slope is too high
00602         //--------------------------------------------------------------------------------------------------
00603   
00604         if(smoothedpoints.size() < 4) return;
00605         DigiMapIter itSmoothedpoints, itSmoothedpointsNext, itSmoothedpointsBegin, itSmoothedpointsEnd; 
00606         itSmoothedpoints=smoothedpoints.begin();
00607         while (smoothedpoints.size() >2 && itSmoothedpoints != --(smoothedpoints.end()) ) { //while we are not at the last point
00608           //if(smoothedpoints.size() <2) break;
00609                 // get info about current and next points
00610                 itSmoothedpointsNext = itSmoothedpoints;
00611                 ++itSmoothedpointsNext;
00612                 float strip1 = itSmoothedpoints->first;
00613                 float strip2 = itSmoothedpointsNext->first;
00614                 float adc1 = itSmoothedpoints->second;
00615                 float adc2 = itSmoothedpointsNext->second;
00616                 float m = (adc2 -adc1)/(strip2 -strip1);
00617         
00618                 if (m>2) { // in case of large positive slope, remove next point and try again from same current point
00619                         smoothedpoints.erase(itSmoothedpointsNext);
00620                 } else if (m<-2) { // in case of large negative slope, remove current point and either...
00621                         // move to next point if we have reached the beginning (post-increment to avoid invalidating pointer during erase) or...
00622                         if(itSmoothedpoints==smoothedpoints.begin()) smoothedpoints.erase(itSmoothedpoints++); 
00623                         // try again from the previous point if we have not reached the beginning
00624                         else smoothedpoints.erase(itSmoothedpoints--); 
00625                 } else { // in case of a flat enough slope, continue on to the next point
00626                         itSmoothedpoints++;
00627                 }
00628                 
00629         }
00630 }
00631 
00632 void inline SiStripAPVRestorer::BaselineFollower(DigiMap& smoothedpoints, std::vector<int16_t>& baseline, const float& median){
00633   
00634   baseline.clear();
00635   DigiMapIter itSmoothedpoints;
00636   
00637   
00638   //if not enough points
00639   if(smoothedpoints.size() < minStripsToFit_){
00640      baseline.insert(baseline.begin(),128, median);
00641   } else {
00642      baseline.insert(baseline.begin(),128, 0);  
00643     
00644     DigiMapIter itSmoothedpointsBegin, itSmoothedpointsEnd;
00645     itSmoothedpointsBegin = smoothedpoints.begin();
00646     itSmoothedpointsEnd = --(smoothedpoints.end());
00647     
00648                                 
00649     uint16_t firstStripFlat = itSmoothedpointsBegin->first;
00650     uint16_t lastStripFlat = itSmoothedpointsEnd->first;
00651     int16_t firstStripFlatADC= itSmoothedpointsBegin->second;
00652     int16_t lastStripFlatADC= itSmoothedpointsEnd->second;
00653    
00654     //adding here the costant line at the extremities 
00655     baseline.erase(baseline.begin(), baseline.begin()+firstStripFlat);
00656     baseline.insert(baseline.begin(), firstStripFlat, firstStripFlatADC);
00657     
00658     baseline.erase(baseline.begin()+lastStripFlat, baseline.end());
00659     baseline.insert(baseline.end(), 128 - lastStripFlat, lastStripFlatADC);
00660     
00661     
00662     //IMPORTANT: the itSmoothedpointsEnd should be at least smaller than smoothedpoints.end() -1
00663     for(itSmoothedpoints = itSmoothedpointsBegin; itSmoothedpoints != itSmoothedpointsEnd; ++itSmoothedpoints){  
00664       DigiMapIter itSmoothedpointsNext = itSmoothedpoints;
00665       ++itSmoothedpointsNext;
00666       float strip1 = itSmoothedpoints->first;
00667       float strip2 = itSmoothedpointsNext->first;
00668       float adc1 = itSmoothedpoints->second;
00669       float adc2 = itSmoothedpointsNext->second;
00670      
00671       baseline[strip1] = adc1;
00672       baseline[strip2] = adc2;
00673       float m = (adc2 -adc1)/(strip2 -strip1);
00674       uint16_t itStrip = strip1 +1;
00675       float stripadc = adc1 + m; 
00676       while(itStrip < strip2){
00677                 baseline[itStrip] = stripadc;
00678                 ++itStrip;
00679                 stripadc+=m;
00680       }
00681       
00682     }
00683     
00684   }
00685 }
00686 
00687 
00688 bool SiStripAPVRestorer::CheckBaseline(const std::vector<int16_t> &baseline) const
00689 {
00690         // The Savitzky-Golay (S-G) filter of any left length nL, right
00691         // length nR, and order m, and with an optional opt equals the
00692         // derivative order (0 for the magnitude, 1 for the first
00693         // derivative, etc.) can be calculated using the following
00694         // Mathematica code:
00695         //
00696         // SavitzkyGolay[m_?IntegerQ, {nL_?IntegerQ, nR_?IntegerQ},
00697         //   opt___] := Module[
00698         //   {a, d},
00699         //   d = If[opt === Null, 0, If[IntegerQ[opt], opt, 0]]; 
00700         //   a = Table[
00701         //     If[i == 0 && j == 0, 1, i^j], {i, -nL, nR}, {j, 0, 
00702         //      m}]; (Inverse[Transpose[a].a].Transpose[a])[[d + 1]]];
00703         //
00704         // The following coefficients can be then calculated by:
00705         //
00706         // N[Join[Table[SavitzkyGolay[2, {k, 16}], {k, 0, 16}], 
00707         //   Table[SavitzkyGolay[2, {16, k}], {k, 15, 0, -1}]]]
00708 
00709         // nLR = max(nL, nR)
00710         static const size_t savitzky_golay_n_l_r = 16;
00711         static const float savitzky_golay_coefficient
00712                 [2 * savitzky_golay_n_l_r + 1][2 * savitzky_golay_n_l_r + 1] = {
00713                 { 0.422085, 0.325077, 0.23839, 0.162023, 0.0959752, 0.0402477,
00714                   -0.00515996, -0.0402477, -0.0650155, -0.0794634, -0.0835913,
00715                   -0.0773994, -0.0608875, -0.0340557, 0.00309598, 0.0505676,
00716                   0.108359 },
00717                 { 0.315789, 0.254902, 0.19969, 0.150155, 0.106295, 0.0681115,
00718                   0.0356037, 0.00877193, -0.0123839, -0.0278638, -0.0376677,
00719                   -0.0417957, -0.0402477, -0.0330237, -0.0201238, -0.00154799,
00720                   0.0227038, 0.0526316 },
00721                 { 0.234586, 0.198496, 0.165207, 0.134719, 0.107032, 0.0821465,
00722                   0.0600619, 0.0407784, 0.024296, 0.0106148, -0.000265369,
00723                   -0.00834439, -0.0136223, -0.0160991, -0.0157747, -0.0126493,
00724                   -0.00672269, 0.00200501, 0.0135338 },
00725                 { 0.172078, 0.153076, 0.135099, 0.118148, 0.102221, 0.0873206,
00726                   0.073445, 0.0605947, 0.0487697, 0.0379699, 0.0281955,
00727                   0.0194463, 0.0117225, 0.00502392, -0.000649351, -0.00529733,
00728                   -0.00892003, -0.0115174, -0.0130895, -0.0136364 },
00729                 { 0.123659, 0.116431, 0.109144, 0.101798, 0.0943921, 0.0869268,
00730                   0.0794021, 0.0718179, 0.0641743, 0.0564712, 0.0487087,
00731                   0.0408868, 0.0330054, 0.0250646, 0.0170644, 0.00900473,
00732                   0.000885613, -0.00729294, -0.0155309, -0.0238283,
00733                   -0.0321852 },
00734                 { 0.0859684, 0.0868154, 0.0869565, 0.0863919, 0.0851214,
00735                   0.0831451, 0.080463, 0.0770751, 0.0729814, 0.0681818,
00736                   0.0626765, 0.0564653, 0.0495483, 0.0419255, 0.0335968,
00737                   0.0245624, 0.0148221, 0.00437606, -0.00677583, -0.0186335,
00738                   -0.0311971, -0.0444664 },
00739                 { 0.0565217, 0.0628458, 0.0680971, 0.0722756, 0.0753811,
00740                   0.0774139, 0.0783738, 0.0782609, 0.0770751, 0.0748165,
00741                   0.071485, 0.0670807, 0.0616036, 0.0550536, 0.0474308,
00742                   0.0387352, 0.0289667, 0.0181254, 0.00621118, -0.00677583,
00743                   -0.0208357, -0.0359684, -0.0521739 },
00744                 { 0.0334615, 0.0434281, 0.0521329, 0.0595759, 0.0657571,
00745                   0.0706765, 0.0743341, 0.07673, 0.0778641, 0.0777364,
00746                   0.0763469, 0.0736957, 0.0697826, 0.0646078, 0.0581712,
00747                   0.0504728, 0.0415126, 0.0312907, 0.0198069, 0.00706142,
00748                   -0.00694588, -0.022215, -0.0387458, -0.0565385 },
00749                 { 0.0153846, 0.0276923, 0.0386622, 0.0482943, 0.0565886,
00750                   0.0635452, 0.0691639, 0.0734448, 0.076388, 0.0779933,
00751                   0.0782609, 0.0771906, 0.0747826, 0.0710368, 0.0659532,
00752                   0.0595318, 0.0517726, 0.0426756, 0.0322408, 0.0204682,
00753                   0.00735786, -0.0070903, -0.0228763, -0.04, -0.0584615 },
00754                 { 0.001221, 0.0149451, 0.027326, 0.0383639, 0.0480586,
00755                   0.0564103, 0.0634188, 0.0690842, 0.0734066, 0.0763858,
00756                   0.078022, 0.078315, 0.077265, 0.0748718, 0.0711355,
00757                   0.0660562, 0.0596337, 0.0518681, 0.0427595, 0.0323077,
00758                   0.0205128, 0.00737485, -0.00710623, -0.0229304, -0.0400977,
00759                   -0.0586081 },
00760                 { -0.00985222, 0.00463138, 0.0178098, 0.029683, 0.0402509,
00761                   0.0495137, 0.0574713, 0.0641236, 0.0694708, 0.0735127,
00762                   0.0762494, 0.0776809, 0.0778073, 0.0766284, 0.0741442,
00763                   0.0703549, 0.0652604, 0.0588607, 0.0511557, 0.0421456,
00764                   0.0318302, 0.0202097, 0.0072839, -0.00694708, -0.0224833,
00765                   -0.0393247, -0.0574713 },
00766                 { -0.0184729, -0.00369458, 0.00984169, 0.0221359, 0.0331881,
00767                   0.0429982, 0.0515662, 0.0588923, 0.0649762, 0.0698181,
00768                   0.073418, 0.0757758, 0.0768915, 0.0767652, 0.0753968,
00769                   0.0727864, 0.0689339, 0.0638394, 0.0575028, 0.0499242,
00770                   0.0411035, 0.0310408, 0.019736, 0.00718917, -0.00659972,
00771                   -0.0216307, -0.0379037, -0.0554187 },
00772                 { -0.025139, -0.0103925, 0.00318873, 0.0156046, 0.0268552,
00773                   0.0369405, 0.0458605, 0.0536151, 0.0602045, 0.0656285,
00774                   0.0698872, 0.0729806, 0.0749086, 0.0756714, 0.0752688,
00775                   0.0737009, 0.0709677, 0.0670692, 0.0620054, 0.0557763,
00776                   0.0483818, 0.039822, 0.0300969, 0.0192065, 0.0071508,
00777                   -0.00607024, -0.0204566, -0.0360083, -0.0527253 },
00778                 { -0.0302419, -0.0157536, -0.00234785, 0.00997537, 0.021216,
00779                   0.0313741, 0.0404497, 0.0484427, 0.0553532, 0.0611811,
00780                   0.0659264, 0.0695892, 0.0721695, 0.0736672, 0.0740823,
00781                   0.0734149, 0.0716649, 0.0688324, 0.0649174, 0.0599198,
00782                   0.0538396, 0.0466769, 0.0384316, 0.0291038, 0.0186934,
00783                   0.00720046, -0.00537502, -0.0190331, -0.0337736,
00784                   -0.0495968 },
00785                 { -0.0340909, -0.0200147, -0.006937, 0.00514208, 0.0162226,
00786                   0.0263045, 0.0353878, 0.0434725, 0.0505587, 0.0566463,
00787                   0.0617353, 0.0658257, 0.0689175, 0.0710107, 0.0721054,
00788                   0.0722014, 0.0712989, 0.0693978, 0.0664981, 0.0625999,
00789                   0.057703, 0.0518076, 0.0449135, 0.0370209, 0.0281297,
00790                   0.01824, 0.0073516, -0.00453534, -0.0174209, -0.031305,
00791                   -0.0461877 },
00792                 { -0.0369318, -0.0233688, -0.0107221, 0.00100806, 0.0118218,
00793                   0.0217192, 0.0307001, 0.0387647, 0.0459128, 0.0521444,
00794                   0.0574597, 0.0618585, 0.0653409, 0.0679069, 0.0695565,
00795                   0.0702896, 0.0701063, 0.0690066, 0.0669905, 0.0640579,
00796                   0.0602089, 0.0554435, 0.0497617, 0.0431635, 0.0356488,
00797                   0.0272177, 0.0178702, 0.0076063, -0.00357405, -0.0156708,
00798                   -0.028684, -0.0426136 },
00799                 { -0.038961, -0.025974, -0.0138249, -0.00251362, 0.00795978,
00800                   0.0175953, 0.026393, 0.0343527, 0.0414747, 0.0477587,
00801                   0.0532049, 0.0578132, 0.0615836, 0.0645161, 0.0666108,
00802                   0.0678676, 0.0682866, 0.0678676, 0.0666108, 0.0645161,
00803                   0.0615836, 0.0578132, 0.0532049, 0.0477587, 0.0414747,
00804                   0.0343527, 0.026393, 0.0175953, 0.00795978, -0.00251362,
00805                   -0.0138249, -0.025974, -0.038961 },
00806                 { -0.0426136, -0.028684, -0.0156708, -0.00357405, 0.0076063,
00807                   0.0178702, 0.0272177, 0.0356488, 0.0431635, 0.0497617,
00808                   0.0554435, 0.0602089, 0.0640579, 0.0669905, 0.0690066,
00809                   0.0701063, 0.0702896, 0.0695565, 0.0679069, 0.0653409,
00810                   0.0618585, 0.0574597, 0.0521444, 0.0459128, 0.0387647,
00811                   0.0307001, 0.0217192, 0.0118218, 0.00100806, -0.0107221,
00812                   -0.0233688, -0.0369318 },
00813                 { -0.0461877, -0.031305, -0.0174209, -0.00453534, 0.0073516,
00814                   0.01824, 0.0281297, 0.0370209, 0.0449135, 0.0518076,
00815                   0.057703, 0.0625999, 0.0664981, 0.0693978, 0.0712989,
00816                   0.0722014, 0.0721054, 0.0710107, 0.0689175, 0.0658257,
00817                   0.0617353, 0.0566463, 0.0505587, 0.0434725, 0.0353878,
00818                   0.0263045, 0.0162226, 0.00514208, -0.006937, -0.0200147,
00819                   -0.0340909 },
00820                 { -0.0495968, -0.0337736, -0.0190331, -0.00537502, 0.00720046,
00821                   0.0186934, 0.0291038, 0.0384316, 0.0466769, 0.0538396,
00822                   0.0599198, 0.0649174, 0.0688324, 0.0716649, 0.0734149,
00823                   0.0740823, 0.0736672, 0.0721695, 0.0695892, 0.0659264,
00824                   0.0611811, 0.0553532, 0.0484427, 0.0404497, 0.0313741,
00825                   0.021216, 0.00997537, -0.00234785, -0.0157536, -0.0302419 },
00826                 { -0.0527253, -0.0360083, -0.0204566, -0.00607024, 0.0071508,
00827                   0.0192065, 0.0300969, 0.039822, 0.0483818, 0.0557763,
00828                   0.0620054, 0.0670692, 0.0709677, 0.0737009, 0.0752688,
00829                   0.0756714, 0.0749086, 0.0729806, 0.0698872, 0.0656285,
00830                   0.0602045, 0.0536151, 0.0458605, 0.0369405, 0.0268552,
00831                   0.0156046, 0.00318873, -0.0103925, -0.025139 },
00832                 { -0.0554187, -0.0379037, -0.0216307, -0.00659972, 0.00718917,
00833                   0.019736, 0.0310408, 0.0411035, 0.0499242, 0.0575028,
00834                   0.0638394, 0.0689339, 0.0727864, 0.0753968, 0.0767652,
00835                   0.0768915, 0.0757758, 0.073418, 0.0698181, 0.0649762,
00836                   0.0588923, 0.0515662, 0.0429982, 0.0331881, 0.0221359,
00837                   0.00984169, -0.00369458, -0.0184729 },
00838                 { -0.0574713, -0.0393247, -0.0224833, -0.00694708, 0.0072839,
00839                   0.0202097, 0.0318302, 0.0421456, 0.0511557, 0.0588607,
00840                   0.0652604, 0.0703549, 0.0741442, 0.0766284, 0.0778073,
00841                   0.0776809, 0.0762494, 0.0735127, 0.0694708, 0.0641236,
00842                   0.0574713, 0.0495137, 0.0402509, 0.029683, 0.0178098,
00843                   0.00463138, -0.00985222 },
00844                 { -0.0586081, -0.0400977, -0.0229304, -0.00710623, 0.00737485,
00845                   0.0205128, 0.0323077, 0.0427595, 0.0518681, 0.0596337,
00846                   0.0660562, 0.0711355, 0.0748718, 0.077265, 0.078315,
00847                   0.078022, 0.0763858, 0.0734066, 0.0690842, 0.0634188,
00848                   0.0564103, 0.0480586, 0.0383639, 0.027326, 0.0149451,
00849                   0.001221 },
00850                 { -0.0584615, -0.04, -0.0228763, -0.0070903, 0.00735786,
00851                   0.0204682, 0.0322408, 0.0426756, 0.0517726, 0.0595318,
00852                   0.0659532, 0.0710368, 0.0747826, 0.0771906, 0.0782609,
00853                   0.0779933, 0.076388, 0.0734448, 0.0691639, 0.0635452,
00854                   0.0565886, 0.0482943, 0.0386622, 0.0276923, 0.0153846 },
00855                 { -0.0565385, -0.0387458, -0.022215, -0.00694588, 0.00706142,
00856                   0.0198069, 0.0312907, 0.0415126, 0.0504728, 0.0581712,
00857                   0.0646078, 0.0697826, 0.0736957, 0.0763469, 0.0777364,
00858                   0.0778641, 0.07673, 0.0743341, 0.0706765, 0.0657571,
00859                   0.0595759, 0.0521329, 0.0434281, 0.0334615 },
00860                 { -0.0521739, -0.0359684, -0.0208357, -0.00677583, 0.00621118,
00861                   0.0181254, 0.0289667, 0.0387352, 0.0474308, 0.0550536,
00862                   0.0616036, 0.0670807, 0.071485, 0.0748165, 0.0770751,
00863                   0.0782609, 0.0783738, 0.0774139, 0.0753811, 0.0722756,
00864                   0.0680971, 0.0628458, 0.0565217 },
00865                 { -0.0444664, -0.0311971, -0.0186335, -0.00677583, 0.00437606,
00866                   0.0148221, 0.0245624, 0.0335968, 0.0419255, 0.0495483,
00867                   0.0564653, 0.0626765, 0.0681818, 0.0729814, 0.0770751,
00868                   0.080463, 0.0831451, 0.0851214, 0.0863919, 0.0869565,
00869                   0.0868154, 0.0859684 },
00870                 { -0.0321852, -0.0238283, -0.0155309, -0.00729294, 0.000885613,
00871                   0.00900473, 0.0170644, 0.0250646, 0.0330054, 0.0408868,
00872                   0.0487087, 0.0564712, 0.0641743, 0.0718179, 0.0794021,
00873                   0.0869268, 0.0943921, 0.101798, 0.109144, 0.116431,
00874                   0.123659 },
00875                 { -0.0136364, -0.0130895, -0.0115174, -0.00892003, -0.00529733,
00876                   -0.000649351, 0.00502392, 0.0117225, 0.0194463, 0.0281955,
00877                   0.0379699, 0.0487697, 0.0605947, 0.073445, 0.0873206,
00878                   0.102221, 0.118148, 0.135099, 0.153076, 0.172078 },
00879                 { 0.0135338, 0.00200501, -0.00672269, -0.0126493, -0.0157747,
00880                   -0.0160991, -0.0136223, -0.00834439, -0.000265369, 0.0106148,
00881                   0.024296, 0.0407784, 0.0600619, 0.0821465, 0.107032,
00882                   0.134719, 0.165207, 0.198496, 0.234586 },
00883                 { 0.0526316, 0.0227038, -0.00154799, -0.0201238, -0.0330237,
00884                   -0.0402477, -0.0417957, -0.0376677, -0.0278638, -0.0123839,
00885                   0.00877193, 0.0356037, 0.0681115, 0.106295, 0.150155,
00886                   0.19969, 0.254902, 0.315789 },
00887                 { 0.108359, 0.0505676, 0.00309598, -0.0340557, -0.0608875,
00888                   -0.0773994, -0.0835913, -0.0794634, -0.0650155, -0.0402477,
00889                   -0.00515996, 0.0402477, 0.0959752, 0.162023, 0.23839,
00890                   0.325077, 0.422085 }
00891         };
00892 
00893         float filtered_baseline[128];
00894         float filtered_baseline_derivative[127];
00895 
00896         // Zero filtered_baseline
00897         memset(filtered_baseline, 0, 128 * sizeof(float));
00898         // Filter the left edge using (nL, nR) = (0, 16) .. (15, 16) S-G
00899         // filters
00900         for (size_t i = 0; i < savitzky_golay_n_l_r; i++) {
00901                 for (size_t j = 0; j < savitzky_golay_n_l_r + 1 + i; j++) {
00902                         filtered_baseline[i] +=
00903                                 savitzky_golay_coefficient[i][j] * baseline[j];
00904                 }
00905         }
00906         // Filter the middle section using the (nL, nR) = (16, 16) S-G
00907         // filter, while taking advantage of the symmetry to save 16
00908         // multiplications.
00909         for (size_t i = savitzky_golay_n_l_r;
00910                  i < 128 - savitzky_golay_n_l_r; i++) {
00911                 filtered_baseline[i] =
00912                         savitzky_golay_coefficient
00913                         [savitzky_golay_n_l_r][savitzky_golay_n_l_r] * baseline[i];
00914                 for (size_t j = 0; j < savitzky_golay_n_l_r; j++) {
00915                         filtered_baseline[i] +=
00916                                 savitzky_golay_coefficient[savitzky_golay_n_l_r][j] *
00917                                 (baseline[i + j - savitzky_golay_n_l_r] +
00918                                  baseline[i - j + savitzky_golay_n_l_r]);
00919                 }
00920 #if 0
00921                 // Test that the indexing above is correct
00922                 float test = 0;
00923                 for (size_t j = 0; j < 2 * savitzky_golay_n_l_r + 1; j++) {
00924                         test +=
00925                                 savitzky_golay_coefficient[savitzky_golay_n_l_r][j] *
00926                                 baseline[i + j - savitzky_golay_n_l_r];
00927                 }
00928                 // test == filtered_baseline[i] should hold now
00929 #endif
00930         }
00931         // Filter the right edge using (nL, nR) = (16, 15) .. (16, 0) S-G
00932         // filters
00933         for (size_t i = 128 - savitzky_golay_n_l_r; i < 128; i++) {
00934                 for (size_t j = 0; j < 128 - i + savitzky_golay_n_l_r; j++) {
00935                         filtered_baseline[i] +=
00936                                 savitzky_golay_coefficient
00937                                 [2 * savitzky_golay_n_l_r + i + 1 - 128][j] *
00938                                 baseline[i + j - savitzky_golay_n_l_r];
00939                 }
00940         }
00941         // In lieu of a spearate S-G derivative filter, the finite
00942         // difference is used here (since the output is sufficiently
00943         // smooth).
00944         for (size_t i = 0; i < 127; i++) {
00945                 filtered_baseline_derivative[i] =
00946                         filtered_baseline[i + 1] - filtered_baseline[i];
00947         }
00948 
00949         // Calculate the maximum deviation between filtered and unfiltered
00950         // baseline, plus the sum square of the derivative.
00951 
00952         double filtered_baseline_max = 0;
00953         double filtered_baseline_derivative_sum_square = 0;
00954 
00955         for (size_t i = 0; i < 128; i++) {
00956                 const double d = filtered_baseline[i] - baseline[i];
00957 
00958                 filtered_baseline_max =
00959                         std::max(filtered_baseline_max,
00960                                          static_cast<double>(fabs(d)));
00961         }
00962         for (size_t i = 0; i < 127; i++) {
00963                 filtered_baseline_derivative_sum_square +=
00964                         filtered_baseline_derivative[i] *
00965                         filtered_baseline_derivative[i];
00966         }
00967 
00968 #if 0
00969         std::cerr << __FILE__ << ':' << __LINE__ << ": "
00970                           << filtered_baseline_max << ' '
00971                           << filtered_baseline_derivative_sum_square << std::endl;
00972 #endif
00973 
00974         // Apply the cut
00975         return !(filtered_baseline_max >= filteredBaselineMax_ ||
00976                          filtered_baseline_derivative_sum_square >= filteredBaselineDerivativeSumSquare_);
00977 }
00978 
00979 
00980 
00981 
00982 
00983 //Other methods implementation ==============================================
00984 //==========================================================================
00985 
00986 void SiStripAPVRestorer::LoadMeanCMMap(const edm::Event& iEvent){
00987   if(useRealMeanCM_){  
00988         edm::Handle< edm::DetSetVector<SiStripRawDigi> > input;
00989     iEvent.getByLabel("siStripDigis","VirginRaw", input);
00990    this->CreateCMMapRealPed(*input);
00991   } else {
00992     edm::Handle< edm::DetSetVector<SiStripProcessedRawDigi> > inputCM;
00993     iEvent.getByLabel("MEANAPVCM",inputCM);
00994     this->CreateCMMapCMstored(*inputCM);
00995   }
00996 }
00997 
00998 
00999 void SiStripAPVRestorer::CreateCMMapRealPed(const edm::DetSetVector<SiStripRawDigi>& input){
01000   
01001   MeanCMmap_.erase(MeanCMmap_.begin(), MeanCMmap_.end());
01002         
01003  //std::cout<< "===============================================" << std::endl;
01004  
01005  for ( edm::DetSetVector<SiStripRawDigi>::const_iterator 
01006           rawDigis = input.begin(); rawDigis != input.end(); rawDigis++) {
01007          SiStripPedestals::Range detPedestalRange = pedestalHandle->getRange(rawDigis->id);
01008                  std::vector<float> MeanCMDetSet;
01009                  MeanCMDetSet.clear();
01010                 
01011                 for(uint16_t APV = 0; APV < rawDigis->size()/128; ++APV){
01012                         uint16_t MinPed =0;
01013                         for(uint16_t strip = APV*128; strip< (APV+1)*128; ++strip){
01014                           uint16_t ped =  (uint16_t)pedestalHandle->getPed(strip,detPedestalRange);
01015                          
01016                           if(ped < MinPed) MinPed = ped;
01017                         }
01018                         if(MinPed>128) MinPed=128;
01019                         MeanCMDetSet.push_back(MinPed);
01020                        
01021                 }
01022                 MeanCMmap_.insert(std::pair<uint32_t, std::vector<float> >(rawDigis->id,MeanCMDetSet));
01023                 
01024  }
01025 }
01026 
01027 void SiStripAPVRestorer::CreateCMMapCMstored(const edm::DetSetVector<SiStripProcessedRawDigi>& Input){
01028 
01029   MeanCMmap_.erase(MeanCMmap_.begin(), MeanCMmap_.end());
01030   uint32_t detId;
01031   edm::DetSetVector<SiStripProcessedRawDigi>::const_iterator itInput;
01032   edm::DetSet<SiStripProcessedRawDigi>::const_iterator itCM;
01033   std::vector<float> MeanCMNValue;
01034   
01035   for(itInput = Input.begin(); itInput != Input.end(); ++itInput){
01036     detId = itInput->id;
01037     MeanCMNValue.clear();
01038     for(itCM = itInput->begin(); itCM != itInput->end(); ++itCM) MeanCMNValue.push_back(itCM->adc());                   
01039     MeanCMmap_.insert(std::pair<uint32_t, std::vector<float> >(detId,MeanCMNValue));
01040   }
01041 }
01042 
01043 std::vector<bool>& SiStripAPVRestorer::GetAPVFlags(){
01044     apvFlagsBool_.clear();
01045     for(size_t i =0; i < apvFlags_.size(); ++i){
01046       if(apvFlags_[i] != "" && !apvFlagsBoolOverride_[i]) apvFlagsBool_.push_back(true);
01047       else apvFlagsBool_.push_back(false);
01048     }
01049     return apvFlagsBool_;       
01050 }
01051 
01052 
01053 
01054