CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/DQM/SiStripCommissioningAnalysis/src/ApvTimingAlgorithm.cc

Go to the documentation of this file.
00001 #include "DQM/SiStripCommissioningAnalysis/interface/ApvTimingAlgorithm.h"
00002 #include "CondFormats/SiStripObjects/interface/ApvTimingAnalysis.h"
00003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
00004 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 #include "TProfile.h"
00007 #include "TH1.h"
00008 #include <iostream>
00009 #include <iomanip>
00010 #include <cmath>
00011 
00012 using namespace sistrip;
00013 
00014 // ----------------------------------------------------------------------------
00015 // 
00016 ApvTimingAlgorithm::ApvTimingAlgorithm( const edm::ParameterSet & pset, ApvTimingAnalysis* const anal ) 
00017   : CommissioningAlgorithm(anal),
00018     histo_(0,"")
00019 {;}
00020 
00021 // ----------------------------------------------------------------------------
00022 // 
00023 void ApvTimingAlgorithm::extract( const std::vector<TH1*>& histos ) { 
00024 
00025   if ( !anal() ) {
00026     edm::LogWarning(mlCommissioning_)
00027       << "[ApvTimingAlgorithm::" << __func__ << "]"
00028       << " NULL pointer to Analysis object!";
00029     return; 
00030   }
00031   
00032   // Check number of histograms
00033   if ( histos.size() != 1 ) {
00034     anal()->addErrorCode(sistrip::numberOfHistos_);
00035   }
00036   
00037   // Extract FED key from histo title
00038   if ( !histos.empty() ) { anal()->fedKey( extractFedKey( histos.front() ) ); }
00039   
00040   // Extract histograms
00041   std::vector<TH1*>::const_iterator ihis = histos.begin();
00042   for ( ; ihis != histos.end(); ihis++ ) {
00043     
00044     // Check for NULL pointer
00045     if ( !(*ihis) ) { continue; }
00046     
00047     // Check name
00048     SiStripHistoTitle title( (*ihis)->GetName() );
00049     if ( title.runType() != sistrip::APV_TIMING ) {
00050       anal()->addErrorCode(sistrip::unexpectedTask_);
00051       continue;
00052     }
00053     
00054     // Extract timing histo
00055     histo_.first = *ihis;
00056     histo_.second = (*ihis)->GetName();
00057     
00058   }
00059   
00060 }
00061 
00062 // ----------------------------------------------------------------------------
00063 // 
00064 void ApvTimingAlgorithm::analyse() { 
00065 
00066   if ( !anal() ) {
00067     edm::LogWarning(mlCommissioning_)
00068       << "[ApvTimingAlgorithm::" << __func__ << "]"
00069       << " NULL pointer to base Analysis object!";
00070     return; 
00071   }
00072 
00073   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>( anal() );
00074   ApvTimingAnalysis* anal = dynamic_cast<ApvTimingAnalysis*>( tmp );
00075   if ( !anal ) {
00076     edm::LogWarning(mlCommissioning_)
00077       << "[ApvTimingAlgorithm::" << __func__ << "]"
00078       << " NULL pointer to derived Analysis object!";
00079     return; 
00080   }
00081 
00082   if ( !histo_.first ) {
00083     anal->addErrorCode(sistrip::nullPtr_);
00084     return;
00085   }
00086   
00087   TProfile* histo = dynamic_cast<TProfile*>(histo_.first);
00088   if ( !histo ) {
00089     anal->addErrorCode(sistrip::nullPtr_);
00090     return;
00091   }
00092   
00093   // Transfer histogram contents/errors/stats to containers
00094   uint16_t non_zero = 0;
00095   float max = -1. * sistrip::invalid_;
00096   float min = 1. * sistrip::invalid_;
00097   uint16_t nbins = static_cast<uint16_t>( histo->GetNbinsX() );
00098   std::vector<float> bin_contents; 
00099   std::vector<float> bin_errors;
00100   std::vector<float> bin_entries;
00101   bin_contents.reserve( nbins );
00102   bin_errors.reserve( nbins );
00103   bin_entries.reserve( nbins );
00104   for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) {
00105     bin_contents.push_back( histo->GetBinContent(ibin+1) );
00106     bin_errors.push_back( histo->GetBinError(ibin+1) );
00107     bin_entries.push_back( histo->GetBinEntries(ibin+1) );
00108     if ( bin_entries[ibin] ) { 
00109       if ( bin_contents[ibin] > max ) { max = bin_contents[ibin]; }
00110       if ( bin_contents[ibin] < min ) { min = bin_contents[ibin]; }
00111       non_zero++;
00112     }
00113   }
00114   if ( bin_contents.size() < 100 ) { 
00115     anal->addErrorCode(sistrip::numberOfBins_);
00116     return; 
00117   }
00118   
00119   // Calculate range (max-min) and threshold level (range/2)
00120   float threshold = min + ( max - min ) / 2.;
00121   anal->base_   = min;
00122   anal->peak_   = max;
00123   anal->height_ = max - min;
00124   if ( max - min < ApvTimingAnalysis::tickMarkHeightThreshold_ ) {
00125     anal->addErrorCode(sistrip::smallDataRange_);
00126     return; 
00127   }
00128   
00129   // Associate samples with either "tick mark" or "baseline"
00130   std::vector<float> tick;
00131   std::vector<float> base;
00132   for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) { 
00133     if ( bin_entries[ibin] ) {
00134       if ( bin_contents[ibin] < threshold ) { 
00135         base.push_back( bin_contents[ibin] ); 
00136       } else { 
00137         tick.push_back( bin_contents[ibin] ); 
00138       }
00139     }
00140   }
00141   
00142   // Find median level of tick mark and baseline
00143   float tickmark = 0.;
00144   float baseline = 0.;
00145   sort( tick.begin(), tick.end() );
00146   sort( base.begin(), base.end() );
00147   if ( !tick.empty() ) { tickmark = tick[ tick.size()%2 ? tick.size()/2 : tick.size()/2 ]; }
00148   if ( !base.empty() ) { baseline = base[ base.size()%2 ? base.size()/2 : base.size()/2 ]; }
00149   anal->base_   = baseline;
00150   anal->peak_   = tickmark;
00151   anal->height_ = tickmark - baseline;
00152   if ( tickmark - baseline < ApvTimingAnalysis::tickMarkHeightThreshold_ ) {
00153     anal->addErrorCode(sistrip::smallTickMarkHeight_);
00154     return; 
00155   }
00156   
00157   // Find rms spread in "baseline" samples
00158   float mean = 0.;
00159   float mean2 = 0.;
00160   for ( uint16_t ibin = 0; ibin < base.size(); ibin++ ) {
00161     mean += base[ibin];
00162     mean2 += base[ibin] * base[ibin];
00163   }
00164   if ( !base.empty() ) { 
00165     mean = mean / base.size();
00166     mean2 = mean2 / base.size();
00167   } else { 
00168     mean = 0.; 
00169     mean2 = 0.; 
00170   }
00171   float baseline_rms = sqrt( fabs( mean2 - mean*mean ) ); 
00172 
00173 
00174   // Find rising edges (derivative across two bins > threshold) 
00175   std::map<uint16_t,float> edges;
00176   for ( uint16_t ibin = 1; ibin < nbins-1; ibin++ ) {
00177     if ( bin_entries[ibin+1] && 
00178          bin_entries[ibin-1] ) {
00179       float derivative = bin_contents[ibin+1] - bin_contents[ibin-1];
00180       if ( derivative > 3.*baseline_rms ) { 
00181         edges[ibin] = derivative; 
00182       }
00183     }
00184   }
00185   if ( edges.empty() ) {
00186     anal->addErrorCode(sistrip::noRisingEdges_);
00187     return;
00188   }
00189   
00190 
00191   // Iterate through "edges" map
00192   uint16_t max_derivative_bin = sistrip::invalid_;
00193   float max_derivative = -1.*sistrip::invalid_;
00194 
00195   bool found = false;
00196   std::map<uint16_t,float>::iterator iter = edges.begin();
00197   while ( !found && iter != edges.end() ) {
00198 
00199     // Iterate through 50 subsequent samples
00200     bool valid = true;
00201     for ( uint16_t ii = 0; ii < 50; ii++ ) {
00202       uint16_t bin = iter->first + ii;
00203 
00204       // Calc local derivative 
00205       float temp = 0.;
00206       if ( static_cast<uint32_t>(bin)   < 1 ||
00207            static_cast<uint32_t>(bin+1) >= nbins ) { 
00208         valid = false; //@@ require complete plateau is found within histo
00209         anal->addErrorCode(sistrip::incompletePlateau_);
00210         continue; 
00211       }
00212       temp = bin_contents[bin+1] - bin_contents[bin-1];
00213       
00214       // Store max derivative
00215       if ( temp > max_derivative ) {
00216         max_derivative = temp;
00217         max_derivative_bin = bin;
00218       }
00219       
00220       // Check if samples following edge are all "high"
00221       if ( ii > 10 && ii < 40 && bin_entries[bin] &&
00222            bin_contents[bin] < baseline + 5.*baseline_rms ) { 
00223         valid = false; 
00224       }
00225 
00226     }
00227 
00228     // Break from loop if tick mark found
00229     if ( valid ) { found = true; }
00230 
00231     /*
00232     else {
00233       max_derivative = -1.*sistrip::invalid_;
00234       max_derivative_bin = sistrip::invalid_;
00235       //edges.erase(iter);
00236       anal->addErrorCode(sistrip::rejectedCandidate_);
00237     }
00238     */
00239 
00240     iter++; // next candidate
00241     
00242   }
00243 
00244 
00245 
00246   if ( !found ) { //Try tick mark recovery
00247 
00248     max_derivative_bin = sistrip::invalid_;
00249     max_derivative = -1.*sistrip::invalid_;
00250 
00251     // Find rising edges_r (derivative_r across five bins > threshold) 
00252     std::map<uint16_t,float> edges_r;
00253     for ( uint16_t ibin_r = 1; ibin_r < nbins-1; ibin_r++ ) {
00254       if  (bin_entries[ibin_r+4] &&
00255            bin_entries[ibin_r+3] && 
00256            bin_entries[ibin_r+2] && 
00257            bin_entries[ibin_r+1] && 
00258            bin_entries[ibin_r] && 
00259            bin_entries[ibin_r-1]) {
00260         float derivative_r = bin_contents[ibin_r+1] - bin_contents[ibin_r-1];
00261         float derivative_r1 = bin_contents[ibin_r+1] - bin_contents[ibin_r];
00262         float derivative_r2 = bin_contents[ibin_r+2] - bin_contents[ibin_r+1];
00263         float derivative_r3 = bin_contents[ibin_r+3] - bin_contents[ibin_r+2];
00264         
00265         if ( derivative_r > 3.*baseline_rms && 
00266              derivative_r1 > 1.*baseline_rms && 
00267              derivative_r2 > 1.*baseline_rms && 
00268              derivative_r3 > 1.*baseline_rms ) { 
00269           edges_r[ibin_r] = derivative_r; 
00270         }
00271       }
00272     }
00273     if ( edges_r.empty() ) {
00274       anal->addErrorCode(sistrip::noRisingEdges_);
00275       return;
00276     }
00277     
00278     // Iterate through "edges_r" map
00279     float max_derivative_r = -1.*sistrip::invalid_;
00280     
00281     bool found_r = false;
00282     std::map<uint16_t,float>::iterator iter_r = edges_r.begin();
00283     while ( !found_r && iter_r != edges_r.end() ) {
00284       
00285       // Iterate through 50 subsequent samples
00286       bool valid_r = true;
00287       int lowpointcount_r = 0;
00288       const int lowpointallow_r = 25;  //Number of points allowed to fall below threshhold w/o invalidating tick mark
00289       for ( uint16_t ii_r = 0; ii_r < 50; ii_r++ ) {
00290         uint16_t bin_r = iter_r->first + ii_r;
00291         
00292         // Calc local derivative_r 
00293         float temp_r = 0.;
00294         if ( static_cast<uint32_t>(bin_r)   < 1 ||
00295              static_cast<uint32_t>(bin_r+1) >= nbins ) { 
00296           valid_r = false; //@@ require complete plateau is found_r within histo
00297           anal->addErrorCode(sistrip::incompletePlateau_);
00298           continue; 
00299         }
00300         temp_r = bin_contents[bin_r+1] - bin_contents[bin_r-1];
00301         
00302         // Store max derivative_r
00303         if ( temp_r > max_derivative_r && ii_r < 10) {
00304           max_derivative_r = temp_r;
00305           max_derivative = temp_r;
00306           max_derivative_bin = bin_r;
00307         }
00308         
00309         
00310         // Check if majority of samples following edge are all "high"
00311         if ( ii_r > 10 && ii_r < 40 && bin_entries[bin_r] &&
00312              bin_contents[bin_r] < baseline + 5.*baseline_rms ) {
00313           lowpointcount_r++;
00314           if(lowpointcount_r > lowpointallow_r){
00315             valid_r = false; 
00316           }
00317         }
00318       }
00319       
00320       // Break from loop if recovery tick mark found
00321       if ( valid_r ) {
00322         found_r = true; 
00323         found = true;
00324         anal->addErrorCode(sistrip::tickMarkRecovered_);
00325       }
00326       else {
00327         max_derivative_r = -1.*sistrip::invalid_;
00328         max_derivative = -1.*sistrip::invalid_;
00329         max_derivative_bin = sistrip::invalid_;
00330         //edges_r.erase(iter_r);
00331         anal->addErrorCode(sistrip::rejectedCandidate_);
00332       }
00333       
00334       iter_r++; // next candidate
00335       
00336     }
00337   } //End tick mark recovery
00338   
00339   // Record time monitorable and check tick mark height
00340   if ( max_derivative_bin <= sistrip::valid_ ) {
00341     anal->time_ = max_derivative_bin * 25. / 24.;
00342     if ( anal->height_ < ApvTimingAnalysis::tickMarkHeightThreshold_ ) { 
00343       anal->addErrorCode(sistrip::tickMarkBelowThresh_);      
00344     }
00345   } else {
00346     anal->addErrorCode(sistrip::missingTickMark_);
00347   }
00348   
00349 }