Reference Manual

Go to the documentation of this file.
00001 #include "RecoLocalTracker/SiStripClusterizer/interface/OldThreeThresholdAlgorithm.h"
00002 #include "sstream"
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 
00005 #define PATCH_FOR_DIGIS_DUPLICATION
00006 
00007 void OldThreeThresholdAlgorithm::initialize(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 OldThreeThresholdAlgorithm::clusterizeDetUnit(const edm::DetSet<SiStripDigi>    & digis, edmNew::DetSetVector<SiStripCluster>::FastFiller & output) {
00015   clusterizeDetUnit_(digis,output);
00016 }
00017 void OldThreeThresholdAlgorithm::clusterizeDetUnit(const edmNew::DetSet<SiStripDigi> & digis, edmNew::DetSetVector<SiStripCluster>::FastFiller & output) {
00018   clusterizeDetUnit_(digis,output);
00019 }
00020 
00021 template<typename InputDetSet>
00022 void OldThreeThresholdAlgorithm::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") << "[OldThreeThresholdAlgorithm::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 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00136     bool isDigiListBad=false;
00137     int16_t oldStrip=-1;
00138 #endif
00139 
00140     for (itest=ibeg; itest<=iend; itest++) {
00141       float channelNoise = noiseHandle_->getNoise(itest->strip(),detNoiseRange);
00142       bool IsBadChannel = qualityHandle_->IsStripBad(detQualityRange,itest->strip());
00143 
00144 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00145       if(itest->strip()==oldStrip){
00146         isDigiListBad=true;
00147         printPatchError=true;
00148         countErrors++;
00149         break;
00150       }
00151       oldStrip=itest->strip();
00152 #endif
00153 
00154 
00155 #ifdef DEBUG_SiStripClusterizer_
00156       if(edm::isDebugEnabled())
00157         ss << "\nLooking at cluster digis:\n\t\t detID " << detID << " digis " << itest->strip()  
00158            << " adc " << itest->adc() << " channelNoise " << channelNoise 
00159            <<  " IsBadChannel from SiStripQuality " << IsBadChannel;
00160 #endif
00161       
00162  
00163       //check for consecutive digis
00164       if (itest!=ibeg && itest->strip()-(itest-1)->strip()!=1){
00165         //digi *(itest-1) and *itest are not consecutive: create an equivalent number of Digis with zero amp
00166         for (int j=(itest-1)->strip()+1;j<itest->strip();j++){
00167           cluster_digis_.push_back(SiStripDigi(j,0)); //if strip bad or under threshold set StripDigi.adc_=0  
00168 #ifdef DEBUG_SiStripClusterizer_
00169           ss << "\n\t\tHole added: detID " << detID << " digis " << j << " adc 0 ";
00170 #endif
00171         }
00172       }
00173       if (!IsBadChannel && itest->adc() >= static_cast<int>( channelThresholdInNoiseSigma()*channelNoise)) {
00174 
00175         float gainFactor  = gainHandle_->getStripGain(itest->strip(), detGainRange);
00176 
00177         // Begin of Change done by Loic Quertenmont
00178         // Protection if the Charge/Gain brings the charge to be bigger 
00179         // than 255 (largest integer containable in uint8)
00180         // Also, Gain is applied only if the strip is not saturating.
00181         // Same convention as in SimTracker/SiStripDigitizer/src/SiTrivialDigitalConverter.cc
00182 
00183         float stripCharge=(static_cast<float>(itest->adc()));
00184 
00185         // //dummy DEBUGG
00186         // float stripCharge;
00187         // for (uint16_t myadc=0; myadc<=255; myadc++){
00188 
00189         //   stripCharge=(static_cast<float>(myadc));
00190 
00191         //   gainFactor=0.73;
00192         // //ENDDEBUGG
00193 
00194         //correct for gain only non-truncated channel charge. Throw exception if channel charge exceeding 255 ADC counts is found
00195  
00196         if(stripCharge<254){
00197           stripCharge /= gainFactor;  
00198               
00199           if(stripCharge>511.5){stripCharge=255;}
00200           else if(stripCharge>253.5){stripCharge=254;}
00201         }  
00202         else if(stripCharge>255){
00203           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";
00204 
00205         }
00206 
00207         //  //dummy DEBUGG
00208         //   std::vector<SiStripDigi> b; b.push_back(SiStripDigi(itest->strip(), static_cast<uint8_t>(stripCharge+0.5)));; 
00209         //   SiStripCluster c( detID, SiStripCluster::SiStripDigiRange( b.begin(),b.end()));
00210         //   edm::LogInfo("MYTEST") << "myadc="<<myadc<<" stripCharge="<<stripCharge<<" stored charge="<< (unsigned int)(c.amplitudes()[0])<< std::endl;
00211         //  edm::LogInfo("MYTEST") << "myadc="<< (unsigned int)(static_cast<uint8_t>(255.5));
00212         // } 
00213         //   //ENDDEBUGG
00214 
00215         // End of Change done by Loic Quertenmont
00216 
00217 
00218 
00219         charge += stripCharge;
00220         sigmaNoise2 += channelNoise*channelNoise/(gainFactor*gainFactor);
00221       
00222         cluster_digis_.push_back(SiStripDigi(itest->strip(), static_cast<uint8_t>(stripCharge+0.5)));
00223       } else {
00224         cluster_digis_.push_back(SiStripDigi(itest->strip(),0)); //if strip bad or under threshold set SiStripDigi.adc_=0
00225         
00226 #ifdef DEBUG_SiStripClusterizer_
00227         if(edm::isDebugEnabled())
00228           ss << "\n\t\tBad or under threshold digis: detID " << detID  << " digis " << itest->strip()  
00229              << " adc " << itest->adc() << " channelNoise " << channelNoise 
00230              <<  " IsBadChannel from SiStripQuality " << IsBadChannel;
00231 #endif
00232       }
00233     }
00234     float sigmaNoise = sqrt(sigmaNoise2);
00235     
00236 #ifdef PATCH_FOR_DIGIS_DUPLICATION
00237     if (charge >= clusterThresholdInNoiseSigma()*sigmaNoise && !isDigiListBad) {
00238 #else
00239       if (charge >= clusterThresholdInNoiseSigma()*sigmaNoise ) {
00240 #endif
00241 
00242 #ifdef DEBUG_SiStripClusterizer_
00243         if(edm::isDebugEnabled())
00244           ss << "\n\t\tCluster accepted :)";
00245 #endif
00246         output.push_back( SiStripCluster( detID, SiStripCluster::SiStripDigiRange( cluster_digis_.begin(),
00247                                                                                    cluster_digis_.end())));
00248         if (_setDetId) output.back().setId( detID);
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") << "[OldThreeThresholdAlgorithm::clusterizeDetUnit] \n There are errors in " << countErrors << "  clusters due to not unique digis ";
00261 #endif
00262 
00263 #ifdef DEBUG_SiStripClusterizer_
00264     LogDebug("SiStripClusterizer") << "[OldThreeThresholdAlgorithm::clusterizeDetUnit] \n" << ss.str();
00265 #endif
00266   }
00267 
00268
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoLocalTracker/SiStripClusterizer/src/OldThreeThresholdAlgorithm.cc
