Analysis for timing run using APV tick marks. More...

#include <ApvTimingAlgorithm.h>

Inheritance diagram for ApvTimingAlgorithm:

Public Member Functions
	ApvTimingAlgorithm (const edm::ParameterSet &pset, ApvTimingAnalysis *const )

const Histo &	histo () const

virtual	~ApvTimingAlgorithm ()

Public Member Functions inherited from CommissioningAlgorithm
void	analysis (const std::vector< TH1 * > &)

	CommissioningAlgorithm (CommissioningAnalysis *const )

	CommissioningAlgorithm ()

virtual	~CommissioningAlgorithm ()

Private Member Functions
void	analyse ()

	ApvTimingAlgorithm ()

void	extract (const std::vector< TH1 * > &)

Private Attributes
Histo	histo_

Additional Inherited Members
Public Types inherited from CommissioningAlgorithm
typedef std::pair< TH1 *, std::string >	Histo

Protected Member Functions inherited from CommissioningAlgorithm
CommissioningAnalysis *const	anal () const

uint32_t	extractFedKey (const TH1 *const )

Detailed Description

Analysis for timing run using APV tick marks.

Author: M. Wingham, R.Bainbridge

Definition at line 16 of file ApvTimingAlgorithm.h.

Constructor & Destructor Documentation

ApvTimingAlgorithm::ApvTimingAlgorithm	(	const edm::ParameterSet &	pset,
		ApvTimingAnalysis * const	anal
	)

Definition at line 16 of file ApvTimingAlgorithm.cc.

   : CommissioningAlgorithm(anal),
     histo_(0,"")
 {;}

virtual ApvTimingAlgorithm::~ApvTimingAlgorithm ( )

inlinevirtual

Definition at line 22 of file ApvTimingAlgorithm.h.

References histo().

22 {;}

ApvTimingAlgorithm::ApvTimingAlgorithm ( )

inlineprivate

Private constructor.

Definition at line 30 of file ApvTimingAlgorithm.h.

References analyse(), and extract().

30 {;}

Member Function Documentation

void ApvTimingAlgorithm::analyse ( )

privatevirtual

Performs histogram anaysis.

Implements CommissioningAlgorithm.

Definition at line 64 of file ApvTimingAlgorithm.cc.

References ApvTimingAnalysis::addErrorCode(), CommissioningAlgorithm::anal(), runEdmFileComparison::base, ApvTimingAnalysis::base_, stringResolutionProvider_cfi::bin, funct::derivative(), runEdmFileComparison::found, ApvTimingAnalysis::height_, histo(), histo_, cuy::ii, sistrip::incompletePlateau_, sistrip::invalid_, hpstanc_transforms::max, RecoTauDiscriminantConfiguration::mean, min(), sistrip::missingTickMark_, sistrip::mlCommissioning_, pileupCalc::nbins, sistrip::noRisingEdges_, sistrip::nullPtr_, sistrip::numberOfBins_, ApvTimingAnalysis::peak_, sistrip::rejectedCandidate_, sistrip::smallDataRange_, sistrip::smallTickMarkHeight_, mathSSE::sqrt(), groupFilesInBlocks::temp, electronIdCutBased_cfi::threshold, sistrip::tickMarkBelowThresh_, ApvTimingAnalysis::tickMarkHeightThreshold_, sistrip::tickMarkRecovered_, ApvTimingAnalysis::time_, tmp, and sistrip::valid_.

Referenced by ApvTimingAlgorithm().

                                  { 
 
   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
     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_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_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_);
   }
   
 }

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

privatevirtual

Extracts and organises histograms.

Implements CommissioningAlgorithm.

Definition at line 23 of file ApvTimingAlgorithm.cc.

References CommissioningAnalysis::addErrorCode(), CommissioningAlgorithm::anal(), sistrip::APV_TIMING, CommissioningAlgorithm::extractFedKey(), CommissioningAnalysis::fedKey(), histo_, sistrip::mlCommissioning_, sistrip::numberOfHistos_, fftjetcommon_cfi::title, and sistrip::unexpectedTask_.

Referenced by ApvTimingAlgorithm().

                                                                 { 
 
   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();
     
   }
   
 }