|
|
|
Analysis for timing run using APV tick marks. More...
#include <ApvTimingAlgorithm.h>
Public Member Functions | |
| ApvTimingAlgorithm (const edm::ParameterSet &pset, ApvTimingAnalysis *const ) | |
| const Histo & | histo () const |
| virtual | ~ApvTimingAlgorithm () |
Private Member Functions | |
| void | analyse () |
| ApvTimingAlgorithm () | |
| void | extract (const std::vector< TH1 * > &) |
Private Attributes | |
| Histo | histo_ |
Analysis for timing run using APV tick marks.
Definition at line 16 of file ApvTimingAlgorithm.h.
| ApvTimingAlgorithm::ApvTimingAlgorithm | ( | const edm::ParameterSet & | pset, |
| ApvTimingAnalysis * const | anal | ||
| ) |
Definition at line 16 of file ApvTimingAlgorithm.cc.
: CommissioningAlgorithm(anal), histo_(0,"") {;}
| virtual ApvTimingAlgorithm::~ApvTimingAlgorithm | ( | ) | [inline, virtual] |
Definition at line 22 of file ApvTimingAlgorithm.h.
{;}
| ApvTimingAlgorithm::ApvTimingAlgorithm | ( | ) | [inline, private] |
| void ApvTimingAlgorithm::analyse | ( | ) | [private, virtual] |
Performs histogram anaysis.
Implements CommissioningAlgorithm.
Definition at line 64 of file ApvTimingAlgorithm.cc.
References ApvTimingAnalysis::addErrorCode(), CommissioningAlgorithm::anal(), newFWLiteAna::base, ApvTimingAnalysis::base_, newFWLiteAna::bin, funct::derivative(), prof2calltree::edges, newFWLiteAna::found, ApvTimingAnalysis::height_, histo(), histo_, sistrip::incompletePlateau_, sistrip::invalid_, max(), plotscripts::mean(), min, sistrip::missingTickMark_, sistrip::mlCommissioning_, ExpressReco_HICollisions_FallBack::nbins, sistrip::noRisingEdges_, sistrip::nullPtr_, sistrip::numberOfBins_, ApvTimingAnalysis::peak_, sistrip::rejectedCandidate_, sistrip::smallDataRange_, sistrip::smallTickMarkHeight_, python::multivaluedict::sort(), mathSSE::sqrt(), cond::rpcobtemp::temp, crabWrap::threshold, sistrip::tickMarkBelowThresh_, ApvTimingAnalysis::tickMarkHeightThreshold_, sistrip::tickMarkRecovered_, ApvTimingAnalysis::time_, tmp, TrackValidation_HighPurity_cff::valid, and sistrip::valid_.
{
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_);
}
}
| void ApvTimingAlgorithm::extract | ( | const std::vector< TH1 * > & | histos | ) | [private, virtual] |
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_, indexGen::title, and sistrip::unexpectedTask_.
{
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();
}
}
| const ApvTimingAlgorithm::Histo & ApvTimingAlgorithm::histo | ( | ) | const [inline] |
Container of histogram pointer and title.
Definition at line 47 of file ApvTimingAlgorithm.h.
References histo_.
Referenced by analyse().
{ return histo_; }
Histo ApvTimingAlgorithm::histo_ [private] |
Container of histogram pointer and title.
Definition at line 41 of file ApvTimingAlgorithm.h.
1.7.3