CMSSW Reference Manual

00001 #include "RecoLocalTracker/SiStripClusterizer/interface/ThreeThresholdStripClusterizer.h"
00002 #include "sstream"
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 
00005 #define PATCH_FOR_DIGIS_DUPLICATION
00006 
00007 void ThreeThresholdStripClusterizer::init(const edm::EventSetup& es) {
00008   //Get ESObject 
00009   es.get<SiStripGainRcd>().get(gainHandle_);
00010   es.get<SiStripNoisesRcd>().get(noiseHandle_);
00011   es.get<SiStripQualityRcd>().get(qualityLabel_,qualityHandle_);
00012 }
00013 
00014 void ThreeThresholdStripClusterizer::clusterizeDetUnit(const edm::DetSet<SiStripDigi>    & digis, edmNew::DetSetVector<SiStripCluster>::FastFiller & output) {
00015     clusterizeDetUnit_(digis,output);
00016 }
00017 void ThreeThresholdStripClusterizer::clusterizeDetUnit(const edmNew::DetSet<SiStripDigi> & digis, edmNew::DetSetVector<SiStripCluster>::FastFiller & output) {
00018     clusterizeDetUnit_(digis,output);
00019 }
00020 
00021 template<typename InputDetSet>
00022 void ThreeThresholdStripClusterizer::clusterizeDetUnit_(const InputDetSet& input,
00023                                                         edmNew::DetSetVector<SiStripCluster>::FastFiller& output) {
00024   
00025 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00026   bool printPatchError=false;
00027   int countErrors=0;
00028 #endif
00029   const uint32_t& detID = input.detId();
00030   
00031   if (!qualityHandle_->IsModuleUsable(detID)){
00032 #ifdef DEBUG_SiStripClusterizer_
00033     LogDebug("SiStripClusterizer") << "[ThreeThresholdStripClusterizer::clusterizeDetUnit] detid " << detID << " skipped because flagged NOT Usable in SiStripQuality " << std::endl; 
00034 #endif
00035     return;
00036   }
00037   
00038   typename InputDetSet::const_iterator begin=input.begin();
00039   typename InputDetSet::const_iterator end  =input.end();
00040   
00041   typename InputDetSet::const_iterator ibeg, iend, ihigh, itest;  
00042   ibeg = iend = begin;
00043   cluster_digis_.clear();
00044   cluster_digis_.reserve(10);
00045 
00046   //output.data.reserve( (end - begin)/3 + 1); // FIXME put back in if needed
00047 
00048   SiStripApvGain::Range detGainRange =  gainHandle_->getRange(detID); 
00049   SiStripNoises::Range detNoiseRange = noiseHandle_->getRange(detID);
00050   SiStripQuality::Range detQualityRange = qualityHandle_->getRange(detID);
00051 
00052   //  AboveSeed predicate(seedThresholdInNoiseSigma(),SiStripNoiseService_,detID);
00053   AboveSeed predicate(seedThresholdInNoiseSigma(), noiseHandle_, detNoiseRange, qualityHandle_, detQualityRange);
00054 
00055 #ifdef DEBUG_SiStripClusterizer_
00056   std::stringstream ss;
00057 #endif
00058 
00059   while ( ibeg != end &&
00060           (ihigh = std::find_if( ibeg, end, predicate)) != end) {
00061 
00062 #ifdef DEBUG_SiStripClusterizer_
00063     if(edm::isDebugEnabled())
00064       ss << "\nSeed Channel:\n\t\t detID "<< detID << " digis " << ihigh->strip() 
00065          << " adc " << ihigh->adc() << " is " 
00066          << " channelNoise " << noiseHandle_->getNoise(ihigh->strip(),detNoiseRange) 
00067          <<  " IsBadChannel from SiStripQuality " << qualityHandle_->IsStripBad(detQualityRange,ihigh->strip());
00068 #endif
00069     
00070     // The seed strip is ihigh. Scan up and down from it, finding nearby strips above
00071     // threshold, allowing for some holes. The accepted cluster runs from strip ibeg
00072     // to iend, and itest is the strip under study, not yet accepted.
00073     iend = ihigh;
00074     itest = iend + 1;
00075     while ( itest != end && (itest->strip() - iend->strip() <= max_holes_ + 1 )) {
00076       float channelNoise = noiseHandle_->getNoise(itest->strip(),detNoiseRange);
00077       bool IsBadChannel = qualityHandle_->IsStripBad(detQualityRange,itest->strip());
00078 
00079 #ifdef DEBUG_SiStripClusterizer_
00080       if(edm::isDebugEnabled())
00081         ss << "\nStrips on the right:\n\t\t detID " << detID << " digis " << itest->strip()  
00082            << " adc " << itest->adc() << " " << " channelNoise " << channelNoise 
00083            <<  " IsBadChannel from SiStripQuality " << IsBadChannel;
00084 #endif 
00085       
00086       if (!IsBadChannel && itest->adc() >= static_cast<int>( channelThresholdInNoiseSigma() * channelNoise)) { 
00087         iend = itest;
00088       }
00089       ++itest;
00090     }
00091     //if the next digi after iend is an adiacent bad digi then insert into candidate cluster
00092     itest=iend+1;
00093     if ( itest != end && (itest->strip() - iend->strip() == 1) && qualityHandle_->IsStripBad(detQualityRange,itest->strip()) ){
00094 #ifdef DEBUG_SiStripClusterizer_
00095       if(edm::isDebugEnabled())
00096         ss << "\n\t\tInserted bad strip at the end edge iend->strip()= " << iend->strip() << " itest->strip() = " << itest->strip();
00097 #endif      
00098       
00099       iend++;
00100     }
00101 
00102     ibeg = ihigh;
00103     itest = ibeg - 1;
00104     while ( itest >= begin && (ibeg->strip() - itest->strip() <= max_holes_ + 1 )) {
00105       float channelNoise = noiseHandle_->getNoise(itest->strip(),detNoiseRange);
00106       bool IsBadChannel = qualityHandle_->IsStripBad(detQualityRange,itest->strip());
00107 
00108 #ifdef DEBUG_SiStripClusterizer_
00109       if(edm::isDebugEnabled())
00110         ss << "\nStrips on the left:\n\t\t detID " << detID << " digis " << itest->strip()
00111            << " adc " << itest->adc() << " " << " channelNoise " << channelNoise 
00112            <<  " IsBadChannel from SiStripQuality " << IsBadChannel;
00113 #endif      
00114       
00115       
00116       if (!IsBadChannel && itest->adc() >= static_cast<int>( channelThresholdInNoiseSigma() * channelNoise)) { 
00117         ibeg = itest;
00118       }
00119       --itest;
00120     }
00121     //if the next digi after ibeg is an adiacent bad digi then insert into candidate cluster
00122     itest=ibeg-1;
00123     if ( itest >= begin && (ibeg->strip() - itest->strip() == 1) &&  qualityHandle_->IsStripBad(detQualityRange,itest->strip()) ) {    
00124 #ifdef DEBUG_SiStripClusterizer_
00125       if(edm::isDebugEnabled())
00126         ss << "\nInserted bad strip at the begin edge ibeg->strip()= " << ibeg->strip() << " itest->strip() = " << itest->strip();
00127 #endif      
00128       ibeg--;
00129     }
00130     
00131     float charge = 0;
00132     float sigmaNoise2=0;
00133     //int counts=0;
00134     cluster_digis_.clear();
00135 
00136 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00137     bool isDigiListBad=false;
00138     int16_t oldStrip=-1;
00139 #endif
00140 
00141     for (itest=ibeg; itest<=iend; itest++) {
00142       float channelNoise = noiseHandle_->getNoise(itest->strip(),detNoiseRange);
00143       bool IsBadChannel = qualityHandle_->IsStripBad(detQualityRange,itest->strip());
00144 
00145 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00146       if(itest->strip()==oldStrip){
00147         isDigiListBad=true;
00148         printPatchError=true;
00149         countErrors++;
00150         break;
00151       }
00152       oldStrip=itest->strip();
00153 #endif
00154 
00155 
00156 #ifdef DEBUG_SiStripClusterizer_
00157       if(edm::isDebugEnabled())
00158         ss << "\nLooking at cluster digis:\n\t\t detID " << detID << " digis " << itest->strip()  
00159            << " adc " << itest->adc() << " channelNoise " << channelNoise 
00160            <<  " IsBadChannel from SiStripQuality " << IsBadChannel;
00161 #endif
00162       
00163  
00164       //check for consecutive digis
00165       if (itest!=ibeg && itest->strip()-(itest-1)->strip()!=1){
00166         //digi *(itest-1) and *itest are not consecutive: create an equivalent number of Digis with zero amp
00167         for (int j=(itest-1)->strip()+1;j<itest->strip();j++){
00168           cluster_digis_.push_back(SiStripDigi(j,0)); //if strip bad or under threshold set StripDigi.adc_=0  
00169 #ifdef DEBUG_SiStripClusterizer_
00170           ss << "\n\t\tHole added: detID " << detID << " digis " << j << " adc 0 ";
00171 #endif
00172         }
00173       }
00174       if (!IsBadChannel && itest->adc() >= static_cast<int>( channelThresholdInNoiseSigma()*channelNoise)) {
00175 
00176         float gainFactor  = gainHandle_->getStripGain(itest->strip(), detGainRange);
00177 
00178         // Begin of Change done by Loic Quertenmont
00179         // Protection if the Charge/Gain brings the charge to be bigger 
00180         // than 255 (largest integer containable in uint8)
00181         // Also, Gain is applied only if the strip is not saturating.
00182         // Same convention as in SimTracker/SiStripDigitizer/src/SiTrivialDigitalConverter.cc
00183 
00184         float stripCharge=(static_cast<float>(itest->adc()));
00185 
00186 //      //dummy DEBUGG
00187 //      float stripCharge;
00188 //      for (uint16_t myadc=0; myadc<=255; myadc++){
00189 
00190 //        stripCharge=(static_cast<float>(myadc));
00191 
00192 //        gainFactor=0.73;
00193 //      //ENDDEBUGG
00194 
00195 //correct for gain only non-truncated channel charge. Throw exception if channel charge exceeding 255 ADC counts is found
00196  
00197           if(stripCharge<254){
00198             stripCharge /= gainFactor;    
00199             
00200             if(stripCharge>511.5){stripCharge=255;}
00201             else if(stripCharge>253.5){stripCharge=254;}
00202           }  
00203           else if(stripCharge>255){
00204             throw cms::Exception("LogicError") << "Cluster charge (" << stripCharge << ") out of range. This clustering algorithm should only work with input charges lower than or equal to 255 ADC counts";
00205 
00206           }
00207 
00208 //      //dummy DEBUGG
00209 //        std::vector<SiStripDigi> b; b.push_back(SiStripDigi(itest->strip(), static_cast<uint8_t>(stripCharge+0.5)));; 
00210 //        SiStripCluster c( detID, SiStripCluster::SiStripDigiRange( b.begin(),b.end()));
00211 //        edm::LogInfo("MYTEST") << "myadc="<<myadc<<" stripCharge="<<stripCharge<<" stored charge="<< (unsigned int)(c.amplitudes()[0])<< std::endl;
00212 //        edm::LogInfo("MYTEST") << "myadc="<< (unsigned int)(static_cast<uint8_t>(255.5));
00213 //      } 
00214 //        //ENDDEBUGG
00215 
00216 // End of Change done by Loic Quertenmont
00217 
00218 
00219 
00220         charge += stripCharge;
00221         sigmaNoise2 += channelNoise*channelNoise/(gainFactor*gainFactor);
00222       
00223         cluster_digis_.push_back(SiStripDigi(itest->strip(), static_cast<uint8_t>(stripCharge+0.5)));
00224       } else {
00225         cluster_digis_.push_back(SiStripDigi(itest->strip(),0)); //if strip bad or under threshold set SiStripDigi.adc_=0
00226         
00227 #ifdef DEBUG_SiStripClusterizer_
00228          if(edm::isDebugEnabled())
00229            ss << "\n\t\tBad or under threshold digis: detID " << detID  << " digis " << itest->strip()  
00230               << " adc " << itest->adc() << " channelNoise " << channelNoise 
00231               <<  " IsBadChannel from SiStripQuality " << IsBadChannel;
00232 #endif
00233       }
00234     }
00235     float sigmaNoise = sqrt(sigmaNoise2);
00236     
00237 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00238     if (charge >= clusterThresholdInNoiseSigma()*sigmaNoise && !isDigiListBad) {
00239 #else
00240     if (charge >= clusterThresholdInNoiseSigma()*sigmaNoise ) {
00241 #endif
00242 
00243 #ifdef DEBUG_SiStripClusterizer_
00244        if(edm::isDebugEnabled())
00245          ss << "\n\t\tCluster accepted :)";
00246 #endif
00247       output.push_back( SiStripCluster( detID, SiStripCluster::SiStripDigiRange( cluster_digis_.begin(),
00248                                                                                       cluster_digis_.end())));
00249     } else {
00250 #ifdef DEBUG_SiStripClusterizer_
00251       if(edm::isDebugEnabled())
00252         ss << "\n\t\tCluster rejected :(";
00253 #endif
00254     }
00255     ibeg = iend+1;
00256   }
00257 
00258 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00259   if(printPatchError)
00260     edm::LogError("SiStripClusterizer") << "[ThreeThresholdStripClusterizer::clusterizeDetUnit] \n There are errors in " << countErrors << "  clusters due to not unique digis ";
00261 #endif
00262 
00263 #ifdef DEBUG_SiStripClusterizer_
00264   LogDebug("SiStripClusterizer") << "[ThreeThresholdStripClusterizer::clusterizeDetUnit] \n" << ss.str();
00265 #endif
00266 }
00267 
00268
ThreeThresholdStripClusterizer.cc
