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