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
00033 if ( histos.size() != 1 ) {
00034 anal()->addErrorCode(sistrip::numberOfHistos_);
00035 }
00036
00037
00038 if ( !histos.empty() ) { anal()->fedKey( extractFedKey( histos.front() ) ); }
00039
00040
00041 std::vector<TH1*>::const_iterator ihis = histos.begin();
00042 for ( ; ihis != histos.end(); ihis++ ) {
00043
00044
00045 if ( !(*ihis) ) { continue; }
00046
00047
00048 SiStripHistoTitle title( (*ihis)->GetName() );
00049 if ( title.runType() != sistrip::APV_TIMING ) {
00050 anal()->addErrorCode(sistrip::unexpectedTask_);
00051 continue;
00052 }
00053
00054
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
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
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
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
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
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
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
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
00200 bool valid = true;
00201 for ( uint16_t ii = 0; ii < 50; ii++ ) {
00202 uint16_t bin = iter->first + ii;
00203
00204
00205 float temp = 0.;
00206 if ( static_cast<uint32_t>(bin) < 1 ||
00207 static_cast<uint32_t>(bin+1) >= nbins ) {
00208 valid = false;
00209 anal->addErrorCode(sistrip::incompletePlateau_);
00210 continue;
00211 }
00212 temp = bin_contents[bin+1] - bin_contents[bin-1];
00213
00214
00215 if ( temp > max_derivative ) {
00216 max_derivative = temp;
00217 max_derivative_bin = bin;
00218 }
00219
00220
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
00229 if ( valid ) { found = true; }
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240 iter++;
00241
00242 }
00243
00244
00245
00246 if ( !found ) {
00247
00248 max_derivative_bin = sistrip::invalid_;
00249 max_derivative = -1.*sistrip::invalid_;
00250
00251
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
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
00286 bool valid_r = true;
00287 int lowpointcount_r = 0;
00288 const int lowpointallow_r = 25;
00289 for ( uint16_t ii_r = 0; ii_r < 50; ii_r++ ) {
00290 uint16_t bin_r = iter_r->first + ii_r;
00291
00292
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;
00297 anal->addErrorCode(sistrip::incompletePlateau_);
00298 continue;
00299 }
00300 temp_r = bin_contents[bin_r+1] - bin_contents[bin_r-1];
00301
00302
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
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
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
00331 anal->addErrorCode(sistrip::rejectedCandidate_);
00332 }
00333
00334 iter_r++;
00335
00336 }
00337 }
00338
00339
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 }