Reference Manual

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