CMSSW_4_4_3_patch1/src/DQM/SiStripCommissioningAnalysis/src/ApvTimingAlgorithm.cc

Go to the documentation of this file.

#include "DQM/SiStripCommissioningAnalysis/interface/ApvTimingAlgorithm.h"
#include "CondFormats/SiStripObjects/interface/ApvTimingAnalysis.h"
#include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
#include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "TProfile.h"
#include "TH1.h"
#include <iostream>
#include <iomanip>
#include <cmath>

using namespace sistrip;

// ----------------------------------------------------------------------------
// 
ApvTimingAlgorithm::ApvTimingAlgorithm( const edm::ParameterSet & pset, ApvTimingAnalysis* const anal ) 
  : CommissioningAlgorithm(anal),
    histo_(0,"")
{;}

// ----------------------------------------------------------------------------
// 
void ApvTimingAlgorithm::extract( const std::vector<TH1*>& histos ) { 

  if ( !anal() ) {
    edm::LogWarning(mlCommissioning_)
      << "[ApvTimingAlgorithm::" << __func__ << "]"
      << " NULL pointer to Analysis object!";
    return; 
  }
  
  // Check number of histograms
  if ( histos.size() != 1 ) {
    anal()->addErrorCode(sistrip::numberOfHistos_);
  }
  
  // Extract FED key from histo title
  if ( !histos.empty() ) { anal()->fedKey( extractFedKey( histos.front() ) ); }
  
  // Extract histograms
  std::vector<TH1*>::const_iterator ihis = histos.begin();
  for ( ; ihis != histos.end(); ihis++ ) {
    
    // Check for NULL pointer
    if ( !(*ihis) ) { continue; }
    
    // Check name
    SiStripHistoTitle title( (*ihis)->GetName() );
    if ( title.runType() != sistrip::APV_TIMING ) {
      anal()->addErrorCode(sistrip::unexpectedTask_);
      continue;
    }
    
    // Extract timing histo
    histo_.first = *ihis;
    histo_.second = (*ihis)->GetName();
    
  }
  
}

// ----------------------------------------------------------------------------
// 
void ApvTimingAlgorithm::analyse() { 

  if ( !anal() ) {
    edm::LogWarning(mlCommissioning_)
      << "[ApvTimingAlgorithm::" << __func__ << "]"
      << " NULL pointer to base Analysis object!";
    return; 
  }

  CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>( anal() );
  ApvTimingAnalysis* anal = dynamic_cast<ApvTimingAnalysis*>( tmp );
  if ( !anal ) {
    edm::LogWarning(mlCommissioning_)
      << "[ApvTimingAlgorithm::" << __func__ << "]"
      << " NULL pointer to derived Analysis object!";
    return; 
  }

  if ( !histo_.first ) {
    anal->addErrorCode(sistrip::nullPtr_);
    return;
  }
  
  TProfile* histo = dynamic_cast<TProfile*>(histo_.first);
  if ( !histo ) {
    anal->addErrorCode(sistrip::nullPtr_);
    return;
  }
  
  // Transfer histogram contents/errors/stats to containers
  uint16_t non_zero = 0;
  float max = -1. * sistrip::invalid_;
  float min = 1. * sistrip::invalid_;
  uint16_t nbins = static_cast<uint16_t>( histo->GetNbinsX() );
  std::vector<float> bin_contents; 
  std::vector<float> bin_errors;
  std::vector<float> bin_entries;
  bin_contents.reserve( nbins );
  bin_errors.reserve( nbins );
  bin_entries.reserve( nbins );
  for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) {
    bin_contents.push_back( histo->GetBinContent(ibin+1) );
    bin_errors.push_back( histo->GetBinError(ibin+1) );
    bin_entries.push_back( histo->GetBinEntries(ibin+1) );
    if ( bin_entries[ibin] ) { 
      if ( bin_contents[ibin] > max ) { max = bin_contents[ibin]; }
      if ( bin_contents[ibin] < min ) { min = bin_contents[ibin]; }
      non_zero++;
    }
  }
  if ( bin_contents.size() < 100 ) { 
    anal->addErrorCode(sistrip::numberOfBins_);
    return; 
  }
  
  // Calculate range (max-min) and threshold level (range/2)
  float threshold = min + ( max - min ) / 2.;
  anal->base_   = min;
  anal->peak_   = max;
  anal->height_ = max - min;
  if ( max - min < ApvTimingAnalysis::tickMarkHeightThreshold_ ) {
    anal->addErrorCode(sistrip::smallDataRange_);
    return; 
  }
  
  // Associate samples with either "tick mark" or "baseline"
  std::vector<float> tick;
  std::vector<float> base;
  for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) { 
    if ( bin_entries[ibin] ) {
      if ( bin_contents[ibin] < threshold ) { 
        base.push_back( bin_contents[ibin] ); 
      } else { 
        tick.push_back( bin_contents[ibin] ); 
      }
    }
  }
  
  // Find median level of tick mark and baseline
  float tickmark = 0.;
  float baseline = 0.;
  sort( tick.begin(), tick.end() );
  sort( base.begin(), base.end() );
  if ( !tick.empty() ) { tickmark = tick[ tick.size()%2 ? tick.size()/2 : tick.size()/2 ]; }
  if ( !base.empty() ) { baseline = base[ base.size()%2 ? base.size()/2 : base.size()/2 ]; }
  anal->base_   = baseline;
  anal->peak_   = tickmark;
  anal->height_ = tickmark - baseline;
  if ( tickmark - baseline < ApvTimingAnalysis::tickMarkHeightThreshold_ ) {
    anal->addErrorCode(sistrip::smallTickMarkHeight_);
    return; 
  }
  
  // Find rms spread in "baseline" samples
  float mean = 0.;
  float mean2 = 0.;
  for ( uint16_t ibin = 0; ibin < base.size(); ibin++ ) {
    mean += base[ibin];
    mean2 += base[ibin] * base[ibin];
  }
  if ( !base.empty() ) { 
    mean = mean / base.size();
    mean2 = mean2 / base.size();
  } else { 
    mean = 0.; 
    mean2 = 0.; 
  }
  float baseline_rms = sqrt( fabs( mean2 - mean*mean ) ); 


  // Find rising edges (derivative across two bins > threshold) 
  std::map<uint16_t,float> edges;
  for ( uint16_t ibin = 1; ibin < nbins-1; ibin++ ) {
    if ( bin_entries[ibin+1] && 
         bin_entries[ibin-1] ) {
      float derivative = bin_contents[ibin+1] - bin_contents[ibin-1];
      if ( derivative > 3.*baseline_rms ) { 
        edges[ibin] = derivative; 
      }
    }
  }
  if ( edges.empty() ) {
    anal->addErrorCode(sistrip::noRisingEdges_);
    return;
  }
  

  // Iterate through "edges" map
  uint16_t max_derivative_bin = sistrip::invalid_;
  float max_derivative = -1.*sistrip::invalid_;

  bool found = false;
  std::map<uint16_t,float>::iterator iter = edges.begin();
  while ( !found && iter != edges.end() ) {

    // Iterate through 50 subsequent samples
    bool valid = true;
    for ( uint16_t ii = 0; ii < 50; ii++ ) {
      uint16_t bin = iter->first + ii;

      // Calc local derivative 
      float temp = 0.;
      if ( static_cast<uint32_t>(bin)   < 1 ||
           static_cast<uint32_t>(bin+1) >= nbins ) { 
        valid = false; //@@ require complete plateau is found within histo
        anal->addErrorCode(sistrip::incompletePlateau_);
        continue; 
      }
      temp = bin_contents[bin+1] - bin_contents[bin-1];
      
      // Store max derivative
      if ( temp > max_derivative ) {
        max_derivative = temp;
        max_derivative_bin = bin;
      }
      
      // Check if samples following edge are all "high"
      if ( ii > 10 && ii < 40 && bin_entries[bin] &&
           bin_contents[bin] < baseline + 5.*baseline_rms ) { 
        valid = false; 
      }

    }

    // Break from loop if tick mark found
    if ( valid ) { found = true; }

    /*
    else {
      max_derivative = -1.*sistrip::invalid_;
      max_derivative_bin = sistrip::invalid_;
      //edges.erase(iter);
      anal->addErrorCode(sistrip::rejectedCandidate_);
    }
    */

    iter++; // next candidate
    
  }



  if ( !found ) { //Try tick mark recovery

    max_derivative_bin = sistrip::invalid_;
    max_derivative = -1.*sistrip::invalid_;

    // Find rising edges_r (derivative_r across five bins > threshold) 
    std::map<uint16_t,float> edges_r;
    for ( uint16_t ibin_r = 1; ibin_r < nbins-1; ibin_r++ ) {
      if  (bin_entries[ibin_r+4] &&
           bin_entries[ibin_r+3] && 
           bin_entries[ibin_r+2] && 
           bin_entries[ibin_r+1] && 
           bin_entries[ibin_r] && 
           bin_entries[ibin_r-1]) {
        float derivative_r = bin_contents[ibin_r+1] - bin_contents[ibin_r-1];
        float derivative_r1 = bin_contents[ibin_r+1] - bin_contents[ibin_r];
        float derivative_r2 = bin_contents[ibin_r+2] - bin_contents[ibin_r+1];
        float derivative_r3 = bin_contents[ibin_r+3] - bin_contents[ibin_r+2];
        
        if ( derivative_r > 3.*baseline_rms && 
             derivative_r1 > 1.*baseline_rms && 
             derivative_r2 > 1.*baseline_rms && 
             derivative_r3 > 1.*baseline_rms ) { 
          edges_r[ibin_r] = derivative_r; 
        }
      }
    }
    if ( edges_r.empty() ) {
      anal->addErrorCode(sistrip::noRisingEdges_);
      return;
    }
    
    // Iterate through "edges_r" map
    uint16_t max_derivative_r_bin = sistrip::invalid_;
    float max_derivative_r = -1.*sistrip::invalid_;
    
    bool found_r = false;
    std::map<uint16_t,float>::iterator iter_r = edges_r.begin();
    while ( !found_r && iter_r != edges_r.end() ) {
      
      // Iterate through 50 subsequent samples
      bool valid_r = true;
      int lowpointcount_r = 0;
      const int lowpointallow_r = 25;  //Number of points allowed to fall below threshhold w/o invalidating tick mark
      for ( uint16_t ii_r = 0; ii_r < 50; ii_r++ ) {
        uint16_t bin_r = iter_r->first + ii_r;
        
        // Calc local derivative_r 
        float temp_r = 0.;
        if ( static_cast<uint32_t>(bin_r)   < 1 ||
             static_cast<uint32_t>(bin_r+1) >= nbins ) { 
          valid_r = false; //@@ require complete plateau is found_r within histo
          anal->addErrorCode(sistrip::incompletePlateau_);
          continue; 
        }
        temp_r = bin_contents[bin_r+1] - bin_contents[bin_r-1];
        
        // Store max derivative_r
        if ( temp_r > max_derivative_r && ii_r < 10) {
          max_derivative_r = temp_r;
          max_derivative = temp_r;
          max_derivative_r_bin = bin_r;
          max_derivative_bin = bin_r;
        }
        
        
        // Check if majority of samples following edge are all "high"
        if ( ii_r > 10 && ii_r < 40 && bin_entries[bin_r] &&
             bin_contents[bin_r] < baseline + 5.*baseline_rms ) {
          lowpointcount_r++;
          if(lowpointcount_r > lowpointallow_r){
            valid_r = false; 
          }
        }
      }
      
      // Break from loop if recovery tick mark found
      if ( valid_r ) {
        found_r = true; 
        found = true;
        anal->addErrorCode(sistrip::tickMarkRecovered_);
      }
      else {
        max_derivative_r = -1.*sistrip::invalid_;
        max_derivative = -1.*sistrip::invalid_;
        max_derivative_r_bin = sistrip::invalid_;
        max_derivative_bin = sistrip::invalid_;
        //edges_r.erase(iter_r);
        anal->addErrorCode(sistrip::rejectedCandidate_);
      }
      
      iter_r++; // next candidate
      
    }
  } //End tick mark recovery
  
  // Record time monitorable and check tick mark height
  if ( max_derivative_bin <= sistrip::valid_ ) {
    anal->time_ = max_derivative_bin * 25. / 24.;
    if ( anal->height_ < ApvTimingAnalysis::tickMarkHeightThreshold_ ) { 
      anal->addErrorCode(sistrip::tickMarkBelowThresh_);      
    }
  } else {
    anal->addErrorCode(sistrip::missingTickMark_);
  }
  
}