CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/DQM/SiStripCommissioningAnalysis/src/FedTimingAlgorithm.cc

Go to the documentation of this file.
00001 #include "DQM/SiStripCommissioningAnalysis/interface/FedTimingAlgorithm.h"
00002 #include "CondFormats/SiStripObjects/interface/FedTimingAnalysis.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 FedTimingAlgorithm::FedTimingAlgorithm( const edm::ParameterSet & pset, FedTimingAnalysis* const anal ) 
00017   : CommissioningAlgorithm(anal),
00018     histo_(0,"")
00019 {;}
00020 
00021 // ----------------------------------------------------------------------------
00022 // 
00023 void FedTimingAlgorithm::extract( const std::vector<TH1*>& histos ) { 
00024   
00025   if ( !anal() ) {
00026     edm::LogWarning(mlCommissioning_)
00027       << "[FedTimingAlgorithm::" << __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::FED_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 FedTimingAlgorithm::analyse() { 
00065 
00066   if ( !anal() ) {
00067     edm::LogWarning(mlCommissioning_)
00068       << "[FedTimingAlgorithm::" << __func__ << "]"
00069       << " NULL pointer to base Analysis object!";
00070     return; 
00071   }
00072 
00073   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>( anal() );
00074   FedTimingAnalysis* anal = dynamic_cast<FedTimingAnalysis*>( tmp );
00075   if ( !anal ) {
00076     edm::LogWarning(mlCommissioning_)
00077       << "[FedTimingAlgorithm::" << __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   // Transfer histogram contents/errors/stats to containers
00088   uint16_t non_zero = 0;
00089   float max = -1.e9;
00090   float min =  1.e9;
00091   uint16_t nbins = static_cast<uint16_t>( histo_.first->GetNbinsX() );
00092   std::vector<float> bin_contents; 
00093   std::vector<float> bin_errors;
00094   std::vector<float> bin_entries;
00095   bin_contents.reserve( nbins );
00096   bin_errors.reserve( nbins );
00097   bin_entries.reserve( nbins );
00098   for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) {
00099     bin_contents.push_back( histo_.first->GetBinContent(ibin+1) );
00100     bin_errors.push_back( histo_.first->GetBinError(ibin+1) );
00101     //bin_entries.push_back( histo_.first->GetBinEntries(ibin+1) );
00102     if ( bin_entries[ibin] ) { 
00103       if ( bin_contents[ibin] > max ) { max = bin_contents[ibin]; }
00104       if ( bin_contents[ibin] < min ) { min = bin_contents[ibin]; }
00105       non_zero++;
00106     }
00107   }
00108 
00109   //LogTrace(mlCommissioning_) << " Number of bins with non-zero entries: " << non_zero;
00110   if ( bin_contents.size() < 100 ) { 
00111     anal->addErrorCode(sistrip::numberOfBins_);
00112     return; 
00113   }
00114   
00115   // Calculate range (max-min) and threshold level (range/2)
00116   float range = max - min;
00117   float threshold = min + range / 2.;
00118   if ( range < 50. ) {
00119     anal->addErrorCode(sistrip::smallDataRange_);
00120     return; 
00121   }
00122   //LogTrace(mlCommissioning_) << " ADC samples: max/min/range/threshold: " 
00123   //<< max << "/" << min << "/" << range << "/" << threshold;
00124   
00125   // Associate samples with either "tick mark" or "baseline"
00126   std::vector<float> tick;
00127   std::vector<float> base;
00128   for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) { 
00129     if ( bin_entries[ibin] ) {
00130       if ( bin_contents[ibin] < threshold ) { 
00131         base.push_back( bin_contents[ibin] ); 
00132       } else { 
00133         tick.push_back( bin_contents[ibin] ); 
00134       }
00135     }
00136   }
00137   //LogTrace(mlCommissioning_) << " Number of 'tick mark' samples: " << tick.size() 
00138   //<< " Number of 'baseline' samples: " << base.size();
00139   
00140   // Find median level of tick mark and baseline
00141   float tickmark = 0.;
00142   float baseline = 0.;
00143   sort( tick.begin(), tick.end() );
00144   sort( base.begin(), base.end() );
00145   if ( !tick.empty() ) { tickmark = tick[ tick.size()%2 ? tick.size()/2 : tick.size()/2 ]; }
00146   if ( !base.empty() ) { baseline = base[ base.size()%2 ? base.size()/2 : base.size()/2 ]; }
00147   //LogTrace(mlCommissioning_) << " Tick mark level: " << tickmark << " Baseline level: " << baseline
00148   //<< " Range: " << (tickmark-baseline);
00149   if ( (tickmark-baseline) < 50. ) {
00150     anal->addErrorCode(sistrip::smallDataRange_);
00151     return; 
00152   }
00153   
00154   // Find rms spread in "baseline" samples
00155   float mean = 0.;
00156   float mean2 = 0.;
00157   for ( uint16_t ibin = 0; ibin < base.size(); ibin++ ) {
00158     mean += base[ibin];
00159     mean2 += base[ibin] * base[ibin];
00160   }
00161   if ( !base.empty() ) { 
00162     mean = mean / base.size();
00163     mean2 = mean2 / base.size();
00164   } else { 
00165     mean = 0.; 
00166     mean2 = 0.; 
00167   }
00168   float baseline_rms = 0.;
00169   if (  mean2 > mean*mean ) { baseline_rms = sqrt( mean2 - mean*mean ); }
00170   else { baseline_rms = 0.; }
00171   //LogTrace(mlCommissioning_) << " Spread in baseline samples: " << baseline_rms;
00172   
00173   // Find rising edges (derivative across two bins > range/2) 
00174   std::map<uint16_t,float> edges;
00175   for ( uint16_t ibin = 1; ibin < nbins-1; ibin++ ) {
00176     if ( bin_entries[ibin+1] && 
00177          bin_entries[ibin-1] ) {
00178       float derivative = bin_contents[ibin+1] - bin_contents[ibin-1];
00179       if ( derivative > 5.*baseline_rms ) {
00180         edges[ibin] = derivative;
00181         //LogTrace(mlCommissioning_) << " Found edge #" << edges.size() << " at bin " << ibin 
00182         //<< " and with derivative " << derivative;
00183       }
00184     }
00185   }
00186   
00187   // Iterate through "edges" std::map
00188   bool found = false;
00189   uint16_t deriv_bin = sistrip::invalid_;
00190   float max_deriv = -1.*sistrip::invalid_;
00191   std::map<uint16_t,float>::iterator iter = edges.begin();
00192   while ( !found && iter != edges.end() ) {
00193 
00194     // Iterate through 50 subsequent samples
00195     bool valid = true;
00196     for ( uint16_t ii = 0; ii < 50; ii++ ) {
00197       uint16_t bin = iter->first + ii;
00198 
00199       // Calc local derivative 
00200       float temp_deriv = 0;
00201       if ( static_cast<uint32_t>(bin)   < 1 ||
00202            static_cast<uint32_t>(bin+1) >= nbins ) { continue; }
00203       temp_deriv = bin_contents[bin+1] - bin_contents[bin-1];
00204       
00205       // Store max derivative
00206       if ( temp_deriv > max_deriv ) {
00207         max_deriv = temp_deriv;
00208         deriv_bin = bin;
00209       }
00210 
00211       // Check if samples following edge are all "high"
00212       if ( ii > 10 && ii < 40 && bin_entries[bin] &&
00213            bin_contents[bin] < baseline + 5*baseline_rms ) { valid = false; }
00214 
00215     }
00216 
00217     // Break from loop if tick mark found
00218     if ( valid ) { found = true; }
00219     else {
00220       max_deriv = -1.*sistrip::invalid_;
00221       deriv_bin = sistrip::invalid_;
00222       edges.erase(iter);
00223     }
00224 
00225     iter++;
00226   }
00227   
00228   // Set monitorables (but not PLL coarse and fine here)
00229   if ( !edges.empty() ) {
00230     anal->time_      = deriv_bin;
00231     anal->error_     = 0.;
00232     anal->base_      = baseline;
00233     anal->peak_      = tickmark;
00234     anal->height_    = tickmark - baseline;
00235   } else {
00236     anal->addErrorCode(sistrip::missingTickMark_);
00237     anal->base_   = baseline;
00238     anal->peak_   = tickmark;
00239     anal->height_ = tickmark - baseline;
00240   }
00241   
00242 }