CMSSW Reference Manual

00001 #include "AnalysisDataFormats/SiStripClusterInfo/interface/SiStripClusterInfo.h"
00002 
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004 #include "FWCore/Framework/interface/Event.h"
00005 #include "FWCore/Framework/interface/EventSetup.h"
00006 #include "FWCore/Framework/interface/ESHandle.h"
00007 
00008 #include "DataFormats/Common/interface/Handle.h"
00009 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
00010 #include "DataFormats/SiStripDigi/interface/SiStripRawDigi.h"
00011 #include "DataFormats/SiStripDigi/interface/SiStripDigi.h"
00012 
00013 #include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h"
00014 #include "CondFormats/DataRecord/interface/SiStripPedestalsRcd.h"
00015 
00016 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
00017 
00018 #include "CommonTools/SiStripZeroSuppression/interface/SiStripMedianCommonModeNoiseSubtraction.h"
00019 #include "CommonTools/SiStripZeroSuppression/interface/SiStripTT6CommonModeNoiseSubtraction.h"
00020 
00021 #include <cmath>
00022 
00023 SiStripClusterInfo::SiStripClusterInfo(const SiStripCluster&  cluster, 
00024                                        const edm::EventSetup&      es,
00025                                        std::string CMNSubtractionMode): es_(es) {
00026   SiStripClusterInfo(cluster.geographicalId(),cluster,es,CMNSubtractionMode);
00027 }
00028 
00029 SiStripClusterInfo::SiStripClusterInfo(const uint32_t  cluster_detId,
00030                                        const SiStripCluster&  cluster, 
00031                                        const edm::EventSetup&      es,
00032                                        std::string CMNSubtractionMode):
00033   es_(es){
00034   
00035   cluster_ = & cluster;  
00036   cluster_detId_ = cluster_detId;  
00037 
00038   es.get<SiStripPedestalsRcd>().get(pedestalsHandle_);
00039   es.get<SiStripNoisesRcd>().get(noiseHandle_);
00040   es.get<SiStripGainRcd>().get(gainHandle_);
00041   
00042   
00043   SiStripPedestalsSubtractor_ = new SiStripPedestalsSubtractor();
00044     
00045     
00046     
00047   //------------------------
00048   if ( CMNSubtractionMode== "Median") { 
00049     SiStripCommonModeNoiseSubtractor_ = new SiStripMedianCommonModeNoiseSubtraction();
00050     validCMNSubtraction_ = true;
00051   }
00052   else if ( CMNSubtractionMode== "TT6") { 
00053     //FIXME : SiStripCommonModeNoiseSubtractor_ = new SiStripTT6CommonModeNoiseSubtraction(conf.getParameter<double>("CutToAvoidSignal"));
00054     validCMNSubtraction_ = true;
00055   }
00056   else {
00057     edm::LogError("SiStripClusterInfoProducer") << "[SiStripClusterInfoProducer::SiStripClusterInfoProducer] No valid CommonModeNoiseSubtraction Mode selected, possible CMNSubtractionMode: Median or TT6" << std::endl;
00058     validCMNSubtraction_ = false;
00059   } 
00060     
00061   //------------------------
00062     
00063   
00064   
00065 }
00066 
00067 
00068 SiStripClusterInfo::~SiStripClusterInfo(){
00069   if ( SiStripPedestalsSubtractor_ != NULL) delete SiStripPedestalsSubtractor_;
00070   if (SiStripCommonModeNoiseSubtractor_ != NULL)  delete SiStripCommonModeNoiseSubtractor_;
00071 }
00072 
00073 
00074 float SiStripClusterInfo::getCharge() const {
00075   
00076   float charge_=0;
00077   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();
00078   
00079   for(size_t i=0; i< amplitudes_.size();i++){
00080     if (amplitudes_[i] > 0){  // redundant as always fullfilled for cluster amplitudes
00081       charge_+=amplitudes_[i];
00082     } 
00083   } 
00084   
00085   return charge_ ;
00086   
00087 }
00088 
00089 
00090 uint16_t SiStripClusterInfo::getMaxPosition() const {
00091   return getMaxPosition(cluster_);
00092 }
00093 
00094 uint16_t SiStripClusterInfo::getMaxPosition(const SiStripCluster *cluster) const {
00095   
00096   uint16_t maxPosition_=0;
00097   const std::vector<uint8_t>& amplitudes_ =  cluster->amplitudes();  
00098   float maxCharge_=0;
00099   
00100   for(size_t i=0; i< amplitudes_.size();i++){
00101     if (amplitudes_[i] > 0){ // redundant as always fullfilled for cluster amplitudes
00102       if (maxCharge_<amplitudes_[i]){ 
00103         maxCharge_=amplitudes_[i];
00104         maxPosition_=i;
00105       }
00106     } 
00107   } 
00108   maxPosition_+=cluster->firstStrip();
00109   
00110   return maxPosition_;
00111   
00112 }
00113 
00114 
00115 float SiStripClusterInfo::getMaxCharge() const {
00116   
00117   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();  
00118   float maxCharge_=0;
00119   
00120   for(size_t i=0; i< amplitudes_.size();i++){
00121     if (amplitudes_[i] > 0){
00122       if (maxCharge_<amplitudes_[i]){ 
00123         maxCharge_=amplitudes_[i];
00124       }
00125     } 
00126   } 
00127   return maxCharge_;
00128 }
00129 
00130 
00131 std::vector<float> SiStripClusterInfo::getRawChargeCLR( const edm::DetSet<SiStripRawDigi>&    ds_SiStripRawDigi,
00132                                                               edmNew::DetSet<SiStripCluster>& ds_SiStripCluster,
00133                                                         std::string                           rawDigiLabel ) {
00134   std::vector<float> adcCLR(3*0);
00135   
00136   // Get First Left and First Right
00137   uint32_t neighbourStrips = 1;
00138   std::pair< std::vector<float>,std::vector<float> > rawChargesLR =
00139     this->getRawDigiAmplitudesLR(neighbourStrips, ds_SiStripRawDigi, ds_SiStripCluster, std::string(rawDigiLabel));
00140   
00141   if (edm::isDebugEnabled()){
00142     std::stringstream sssL;
00143     std::stringstream sssR;
00144     // Left
00145     for(std::vector<float>::const_iterator digi_adc_iter=rawChargesLR.first.begin();
00146         digi_adc_iter!=rawChargesLR.first.end();digi_adc_iter++)
00147       sssL << "\n digi adc " << *digi_adc_iter;
00148     LogDebug("SiStripClusterInfo") << " detid " << cluster_detId_  << "\n Max Left RawDigis \n " << sssL.str();
00149     // Right
00150     for(std::vector<float>::const_iterator digi_adc_iter=rawChargesLR.second.begin();
00151         digi_adc_iter!=rawChargesLR.second.end();digi_adc_iter++)
00152       sssR << "\n digi adc " << *digi_adc_iter;
00153     LogDebug("SiStripClusterInfo") << " detid " << cluster_detId_  << "\n Max Right RawDigis \n " << sssR.str();
00154   }
00155   
00156   adcCLR[0] = amplitudeC_;
00157   adcCLR[1] = rawChargesLR.first[0];
00158   adcCLR[2] = rawChargesLR.second[0];
00159   
00160   return adcCLR;
00161 }
00162 
00163 
00164 std::pair<float,float> SiStripClusterInfo::getChargeLR() const {
00165   
00166   float chargeL_=0;
00167   float chargeR_=0; 
00168   uint16_t maxPosition_ = this->getMaxPosition();
00169   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();
00170   
00171   for(size_t i=0; i<amplitudes_.size();i++){
00172     if ((cluster_->firstStrip())+i < maxPosition_) chargeL_+=(float)amplitudes_[i];       
00173     if ((cluster_->firstStrip())+i > maxPosition_) chargeR_+=(float)amplitudes_[i];       
00174   }
00175   
00176   return std::pair<float,float>(chargeL_,chargeR_);
00177   
00178 }
00179 
00180 
00181 std::pair<float,float> SiStripClusterInfo::getChargeLRFirstNeighbour() const {
00182   
00183   uint16_t maxPosition = this->getMaxPosition();
00184   uint16_t firstStrip  = this->getFirstStrip();
00185   const std::vector<uint8_t>& amplitudes =  cluster_->amplitudes();
00186   
00187   float chargeL = ( maxPosition > firstStrip                           ? (float)amplitudes[maxPosition-firstStrip-1] : 0. );
00188   float chargeR = ( maxPosition < firstStrip + (amplitudes.size() - 1) ? (float)amplitudes[maxPosition-firstStrip+1] : 0. );
00189   
00190   //  std::cout << this->getDetId() << ": " << chargeL << " " << chargeR << " width " << amplitudes.size() << std::endl;
00191   
00192   return std::pair<float,float>(chargeL,chargeR);
00193   
00194 }
00195 
00196 
00197 std::vector<float> SiStripClusterInfo::getStripNoises() const {
00198   
00199   std::vector<float>   stripNoises_;
00200   SiStripNoises::Range detNoiseRange = noiseHandle_->getRange(cluster_detId_);  
00201   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();
00202   
00203   for(size_t i=0; i<amplitudes_.size();i++){
00204     
00205     float noise_=noiseHandle_->getNoise(cluster_->firstStrip()+i,detNoiseRange);        
00206     stripNoises_.push_back(noise_);
00207     
00208   } 
00209   
00210   return stripNoises_;
00211   
00212 }
00213 
00214 
00215 std::vector<float> SiStripClusterInfo::getStripNoisesRescaledByGain() const {
00216   
00217   std::vector<float>   stripNoises_;
00218   SiStripNoises::Range detNoiseRange = noiseHandle_->getRange(cluster_detId_);  
00219   SiStripApvGain::Range detGainRange = gainHandle_->getRange(cluster_detId_);   
00220   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();
00221   
00222   for(size_t i=0; i<amplitudes_.size();i++){
00223     
00224     float gain_=gainHandle_->getStripGain(cluster_->firstStrip()+i,detGainRange);
00225     float noise_=noiseHandle_->getNoise(cluster_->firstStrip()+i,detNoiseRange);
00226     noise_=noise_/gain_;
00227     
00228     stripNoises_.push_back(noise_);
00229     
00230   } 
00231   
00232   return stripNoises_;
00233   
00234 }
00235 
00236 
00237 float SiStripClusterInfo::getNoise() {
00238     
00239   int numberOfPosAmplitudes_=0;
00240   float        clusterNoise_=0;  
00241   float       clusterNoise2_=0;
00242   
00243   SiStripNoises::Range detNoiseRange = noiseHandle_->getRange(cluster_detId_);
00244   
00245   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();
00246   
00247   for(size_t i=0; i<amplitudes_.size();i++){
00248     
00249     float noise_=noiseHandle_->getNoise(cluster_->firstStrip()+i,detNoiseRange);
00250     
00251     if (amplitudes_[i]>0){
00252       clusterNoise2_+=noise_*noise_;
00253       numberOfPosAmplitudes_++;  
00254     } 
00255   }   
00256   clusterNoise_= sqrt(clusterNoise2_/numberOfPosAmplitudes_);
00257   
00258   return clusterNoise_;
00259   
00260 }
00261 
00262 
00263 float SiStripClusterInfo::getNoiseRescaledByGain() {
00264     
00265   int numberOfPosAmplitudes_=0;
00266   float       clusterNoise_ =0;  
00267   float       clusterNoise2_=0;  
00268   
00269   SiStripNoises::Range detNoiseRange = noiseHandle_->getRange(cluster_detId_);
00270   SiStripApvGain::Range detGainRange = gainHandle_->getRange(cluster_detId_);   
00271   
00272   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();
00273   
00274   for(size_t i=0; i<amplitudes_.size();i++){
00275     
00276     float gain_=gainHandle_->getStripGain(cluster_->firstStrip()+i,detGainRange);
00277     float noise_=noiseHandle_->getNoise(cluster_->firstStrip()+i,detNoiseRange);
00278     noise_=noise_/gain_;
00279     
00280     if (amplitudes_[i]>0){
00281       clusterNoise2_+=noise_*noise_;
00282       numberOfPosAmplitudes_++;  
00283     } 
00284   }   
00285   clusterNoise_= sqrt(clusterNoise2_/numberOfPosAmplitudes_);
00286   
00287   return clusterNoise_;
00288   
00289 }
00290 
00291 
00292 float SiStripClusterInfo::getNoiseForStripNb(uint16_t istrip) const{
00293   
00294   std::vector<float>   stripNoises_;    
00295   stripNoises_ = this->getStripNoises();
00296   
00297   short strip_of_istrip     = (int) ((cluster_->firstStrip()+istrip)/128);  
00298   short strip_of_firstStrip_ = (int) (cluster_->firstStrip()/128);  
00299   
00300   return stripNoises_[strip_of_istrip-strip_of_firstStrip_];
00301   
00302 }
00303 
00304 
00305 std::vector<float> SiStripClusterInfo::getApvGains() const {
00306   
00307   std::vector<float>   apvGains_;
00308   SiStripApvGain::Range detGainRange = gainHandle_->getRange(cluster_detId_);   
00309   
00310   const std::vector<uint8_t>& amplitudes_ =  cluster_->amplitudes();
00311   
00312   for(size_t i=0; i<amplitudes_.size();i++){    
00313     float gain_=gainHandle_->getStripGain(cluster_->firstStrip()+i,detGainRange);
00314     if (apvGains_.empty())           apvGains_.push_back(gain_);
00315     else if (apvGains_.back()!=gain_) apvGains_.push_back(gain_);         
00316   } 
00317   
00318   return apvGains_;
00319   
00320 }
00321 
00322 
00323 float SiStripClusterInfo::getGainForStripNb(uint16_t istrip) const{
00324   
00325   std::vector<float>   apvGains_;       
00326   apvGains_ = this->getApvGains();
00327   
00328   short apv_of_istrip = (int) ((cluster_->firstStrip()+istrip)/128);  
00329   short apv_of_firstStrip_ = (int) (cluster_->firstStrip()/128);  
00330   
00331   return apvGains_[apv_of_istrip-apv_of_firstStrip_];
00332   
00333 }
00334 
00335 
00336 float SiStripClusterInfo::getSignalOverNoise() {
00337   
00338   return (this->getCharge())/(this->getNoise());
00339   
00340 }
00341 
00342 
00343 float SiStripClusterInfo::getSignalOverNoiseRescaledByGain() {
00344   
00345   return (this->getCharge())/(this->getNoiseRescaledByGain());
00346   
00347 }
00348 
00349 
00350 std::pair< std::vector<float>,std::vector<float> > SiStripClusterInfo::getRawDigiAmplitudesLR(      uint32_t&                          neighbourStripNr, 
00351                                                                                               const edm::DetSetVector<SiStripRawDigi>&    rawDigis_dsv_, 
00352                                                                                                     edmNew::DetSetVector<SiStripCluster>& clusters_dsv_,
00353                                                                                                     std::string                            rawDigiLabel ) {
00354   neighbourStripNr_ = neighbourStripNr;
00355   
00356   const edm::DetSet<SiStripRawDigi> rawDigis_ds_ = rawDigis_dsv_[cluster_detId_];
00357   edmNew::DetSet<SiStripCluster>    clusters_ds_ = clusters_dsv_[cluster_detId_];
00358  
00359   if (rawDigis_dsv_.find(cluster_detId_)!=rawDigis_dsv_.end()){
00360   
00361     rawdigi_algorithm(rawDigis_ds_,clusters_ds_,rawDigiLabel);
00362   
00363     return std::pair< std::vector<float>,std::vector<float> > (amplitudesL_, amplitudesR_);
00364   }    
00365   else {
00366     throw cms::Exception("CorruptedData")
00367       << "[SiStripClusterInfo::getRawDigiAmplitudesLR] reached already end " << std::endl;
00368   }
00369 }
00370 
00371 
00372 std::pair< std::vector<float>,std::vector<float> > SiStripClusterInfo::getRawDigiAmplitudesLR(      uint32_t&                   neighbourStripNr, 
00373                                                                                               const edm::DetSet<SiStripRawDigi>&    rawDigis_ds_, 
00374                                                                                                     edmNew::DetSet<SiStripCluster>& clusters_ds_,
00375                                                                                                     std::string                     rawDigiLabel ) {
00376   neighbourStripNr_ = neighbourStripNr;
00377   
00378   //QUESTION: here no test?
00379   rawdigi_algorithm(rawDigis_ds_,clusters_ds_,rawDigiLabel);
00380   
00381   return std::pair< std::vector<float>,std::vector<float> > (amplitudesL_, amplitudesR_);
00382 }
00383 
00384 
00385 
00386 std::pair< std::vector<float>,std::vector<float> > SiStripClusterInfo::getDigiAmplitudesLR(      uint32_t&                          neighbourStripNr,
00387                                                                                            const edm::DetSetVector<SiStripDigi>&          digis_dsv_,
00388                                                                                                  edmNew::DetSetVector<SiStripCluster>& clusters_dsv_ ){
00389   
00390   neighbourStripNr_ = neighbourStripNr;
00391   
00392   const edm::DetSet<SiStripDigi>    digis_ds_ =   digis_dsv_[cluster_detId_];
00393   edmNew::DetSet<SiStripCluster> clusters_ds_ = clusters_dsv_[cluster_detId_];
00394 
00395   if (digis_dsv_.find(cluster_detId_)!=digis_dsv_.end()){
00396   
00397     digi_algorithm(digis_ds_,clusters_ds_);          
00398     
00399     return std::pair< std::vector<float>,std::vector<float> > (amplitudesL_, amplitudesR_);
00400   }
00401   else {
00402     throw cms::Exception("CorruptedData")
00403       << "[SiStripClusterInfo::getDigiAmplitudesLR] reached already end " << std::endl;
00404   }
00405  
00406 }
00407 
00408 std::pair< std::vector<float>,std::vector<float> > SiStripClusterInfo::getDigiAmplitudesLR(      uint32_t&                  neighbourStripNr, 
00409                                                                                            const edm::DetSet<SiStripDigi>&          digis_ds_, 
00410                                                                                                  edmNew::DetSet<SiStripCluster>& clusters_ds_ ) {
00411   neighbourStripNr_ = neighbourStripNr;
00412   
00413   digi_algorithm(digis_ds_,clusters_ds_);
00414   
00415   return std::pair< std::vector<float>,std::vector<float> > (amplitudesL_, amplitudesR_);
00416 }
00417 
00418 
00419 void SiStripClusterInfo::rawdigi_algorithm(const edm::DetSet<SiStripRawDigi>     rawDigis_ds_,
00420                                                  edmNew::DetSet<SiStripCluster>  clusters_ds_,
00421                                                  std::string                     rawDigiLabel ){
00422   
00423   std::vector<int16_t> vssRd(rawDigis_ds_.size());
00424   
00425   if ( rawDigiLabel == "ProcessedRaw"){
00426     
00427     for(edm::DetSet<SiStripRawDigi>::const_iterator digis_iter=rawDigis_ds_.data.begin();
00428                                                    digis_iter!=rawDigis_ds_.data.end();digis_iter++){
00429       vssRd.push_back(digis_iter->adc());
00430       
00431       if (edm::isDebugEnabled()){
00432         std::stringstream sss;
00433         int idig=0;
00434         
00435         for(std::vector<int16_t>::const_iterator digi_adc_iter=vssRd.begin();
00436                                                 digi_adc_iter!=vssRd.end();digi_adc_iter++)
00437           sss << "\n digi strip " << idig++ << " digi adc " << *digi_adc_iter;
00438         LogDebug("SiStripClusterInfo") << " detid " << cluster_detId_  << " Pedestal subtracted digis \n" << sss.str(); 
00439       }     
00440     }
00441   } else if ( rawDigiLabel == "VirginRaw" ) {
00442     
00443     if (edm::isDebugEnabled()){
00444       std::stringstream sss;
00445       int idig=0;
00446       
00447       for(edm::DetSet<SiStripRawDigi>::const_iterator digis_iter=rawDigis_ds_.data.begin();
00448                                                      digis_iter!=rawDigis_ds_.data.end();digis_iter++)
00449         sss << "\n digi strip " << idig++ << " digi adc " << digis_iter->adc();
00450       LogDebug("SiStripClusterInfo") << " detid " << cluster_detId_  << "\n RawDigis \n" << sss.str();
00451     }
00452     
00453     //Subtract Pedestals
00454     SiStripPedestalsSubtractor_->init(es_);
00455     SiStripPedestalsSubtractor_->subtract(rawDigis_ds_,vssRd);
00456     
00457     if (edm::isDebugEnabled()){
00458       std::stringstream sss;
00459       int idig=0;
00460       
00461       for(std::vector<int16_t>::const_iterator digi_adc_iter=vssRd.begin();
00462                                               digi_adc_iter!=vssRd.end();digi_adc_iter++)
00463         sss << "\n digi strip " << idig++ << " digi adc " << *digi_adc_iter;
00464       LogDebug("SiStripClusterInfo") << " detid " << cluster_detId_  << "\n Pedestal subtracted digis \n " << sss.str();        
00465     }
00466     
00467     //Subtract CMN
00468     if (validCMNSubtraction_){
00469       SiStripCommonModeNoiseSubtractor_->init(es_);
00470       SiStripCommonModeNoiseSubtractor_->subtract(rawDigis_ds_.id,vssRd);
00471       
00472       if (edm::isDebugEnabled()){
00473         std::stringstream sss;
00474         int idig=0;
00475         
00476         for(std::vector<int16_t>::const_iterator digi_adc_iter=vssRd.begin();
00477                                                 digi_adc_iter!=vssRd.end();digi_adc_iter++)
00478           sss << "\n digi strip " << idig++ << " digi adc " << *digi_adc_iter;
00479         LogDebug("SiStripClusterInfo") << " detid " << cluster_detId_  << "\n CMN subtracted digis \n " << sss.str();   
00480       }
00481       
00482     }else{
00483       throw cms::Exception("") 
00484         << "[" << __PRETTY_FUNCTION__<< "] No valid CommonModeNoiseSubtraction Mode selected, possible CMNSubtractionMode: Median or TT6"
00485         << std::endl;
00486     }
00487   } else {
00488     return;
00489   }  
00490   findNeigh("raw",clusters_ds_,vssRd,vssRd);    
00491 }
00492 
00493 
00494 void SiStripClusterInfo::digi_algorithm(const edm::DetSet<SiStripDigi>          digis_ds_,
00495                                               edmNew::DetSet<SiStripCluster> clusters_ds_ ){
00496   
00497   std::vector<int16_t> vstrip;
00498   std::vector<int16_t> vadc;
00499   
00500   //Get list of digis for the current DetId
00501   vstrip.clear();
00502   vadc.clear();
00503   
00504   for(edm::DetSet<SiStripDigi>::const_iterator digis_iter=digis_ds_.data.begin();
00505                                                digis_iter!=digis_ds_.data.end();digis_iter++){
00506     vstrip.push_back(digis_iter->strip());
00507     vadc.push_back(digis_iter->adc());
00508   }
00509   findNeigh("digi",clusters_ds_,vadc,vstrip);
00510 }  
00511 
00512 
00513 
00514 void SiStripClusterInfo::findNeigh(char*                                  mode,
00515                                    edmNew::DetSet<SiStripCluster> clusters_ds_,
00516                                    std::vector<int16_t>&                  vadc,
00517                                    std::vector<int16_t>&                vstrip ){
00518   
00519   // clean the private members
00520   //  before filling with the neighbours
00521   amplitudeC_ = 0;
00522   amplitudesL_.clear();
00523   amplitudesR_.clear();
00524   //
00525   
00526   // Find Digi adiacent to the clusters of this detid
00527   int16_t lastStrip_previousCluster=-1;
00528   int16_t firstStrip_nextCluster=10000;
00529   
00530   for (edmNew::DetSet<SiStripCluster>::const_iterator clusters_iter=clusters_ds_.begin(); 
00531        clusters_iter!=clusters_ds_.end(); clusters_iter++){     
00532     
00533     // Avoid overlapping with neighbour clusters
00534     if (clusters_iter!=clusters_ds_.begin())
00535       lastStrip_previousCluster=(clusters_iter-1)->firstStrip()+(clusters_iter-1)->amplitudes().size() -1;        
00536     if (clusters_iter!=clusters_ds_.end()-1)
00537       firstStrip_nextCluster=(clusters_iter+1)->firstStrip();
00538     
00539     
00540     // Get Gain Range
00541     SiStripApvGain::Range detGainRange = gainHandle_->getRange(clusters_ds_.id());
00542     
00543     // Left or Right with respect to the Central Strip (Cluster Seed)
00544     //    int16_t firstStrip=clusters_iter->firstStrip();
00545     //    int16_t lastStrip=firstStrip + clusters_iter->amplitudes().size() -1;
00546     int16_t firstStrip=getMaxPosition(clusters_iter);
00547     int16_t lastStrip=firstStrip;
00548     //    std::cout << "firstStrip = " << firstStrip << " lastStrip = " << lastStrip << std::endl;
00549     std::vector<int16_t>::iterator ptr;
00550     if (mode=="digi"){
00551       ptr=std::find(vstrip.begin(),vstrip.end(),firstStrip); 
00552       if (ptr==vstrip.end())
00553         throw cms::Exception("") << "\n Expected Digi not found in detid " << clusters_ds_.id() << " strip " << firstStrip << std::endl;
00554     }
00555     else{ 
00556       ptr=vstrip.begin()+firstStrip; // For raw mode vstrip==vadc==vector of digis for all strips in the det
00557     }
00558     
00559     // Central Digi from the ptr to the central strip
00560     int centralPos = ptr-vstrip.begin();
00561     if (mode=="digi")
00562       centralPos=*(vstrip.begin()+centralPos);
00563     //    std::cout << "centralPos = " << centralPos << std::endl;
00564     float gain=gainHandle_->getStripGain(centralPos,detGainRange);
00565     amplitudeC_ = ( (*(vadc.begin()+(ptr-vstrip.begin()))) / gain );
00566     
00567     // Looking at digis before firstStrip         
00568     for (uint16_t istrip=1;istrip<neighbourStripNr_+1;istrip++){
00569       if (istrip>ptr-vstrip.begin()) //avoid underflow
00570         {break;}
00571       if (mode=="digi")
00572         if (firstStrip-istrip!=*(ptr-istrip)) //avoid not contiguous digis
00573           {break;}
00574       if (firstStrip-istrip==lastStrip_previousCluster) //avoid clusters overlapping 
00575         {break;}
00576       int stripPos=ptr-vstrip.begin()-istrip;
00577       if (mode=="digi")
00578         stripPos=*(vstrip.begin()+stripPos);
00579       //      std::cout << "stripPos = " << stripPos << std::endl;
00580       float gain=gainHandle_->getStripGain(stripPos,detGainRange);
00581       amplitudesL_.push_back( (*(vadc.begin()+(ptr-vstrip.begin())-istrip)) / gain );
00582     }
00583     
00584     ptr+=lastStrip-firstStrip;
00585     
00586     // Looking at digis after LastStrip
00587     for (uint16_t istrip=1;istrip<neighbourStripNr_+1;istrip++){
00588       if (istrip>vstrip.end()-ptr-1) //avoid overflow
00589         {break;}
00590       if (mode=="digi")
00591         if (lastStrip+istrip!=*(ptr+istrip)) //avoid not contiguous digis
00592           {break;}
00593       if (lastStrip+istrip==firstStrip_nextCluster) //avoid clusters overlapping 
00594         {break;}
00595       int stripPos=ptr-vstrip.begin()+istrip;
00596       if (mode=="digi")
00597         stripPos=*(vstrip.begin()+stripPos);
00598       //      std::cout << "stripPos = " << stripPos << std::endl;
00599       float gain=gainHandle_->getStripGain(stripPos,detGainRange);
00600       amplitudesR_.push_back( (*(vadc.begin()+(ptr-vstrip.begin())+istrip)) / gain );
00601     }   
00602   }
00603 }
00604 
00605 
00606 
00607 
00608
SiStripClusterInfo.cc
