SiStripCommissioningSources/src/OptoScanTask.cc

0001 #include "DQM/SiStripCommissioningSources/interface/OptoScanTask.h"
0002 #include "DataFormats/SiStripCommon/interface/SiStripConstants.h"
0003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0004 #include "DQMServices/Core/interface/DQMStore.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include <algorithm>
0007 #include <cmath>
0008
0009 #include "DataFormats/SiStripCommon/interface/SiStripFedKey.h"
0010 #include <iomanip>
0011
0012 using namespace sistrip;
0013
0014 // -----------------------------------------------------------------------------
0015 //
0016 OptoScanTask::OptoScanTask(DQMStore* dqm, const FedChannelConnection& conn)
0017     : CommissioningTask(dqm, conn, "OptoScanTask"), opto_() {}
0018
0019 // -----------------------------------------------------------------------------
0020 //
0021 OptoScanTask::~OptoScanTask() {}
0022
0023 // -----------------------------------------------------------------------------
0024 //
0025 void OptoScanTask::book() {
0026   uint16_t nbins = 50;
0027   uint16_t gains = 4;
0028
0029   // Resize "histo sets"
0030   opto_.resize(gains);
0031   for (uint16_t igain = 0; igain < opto_.size(); igain++) {
0032     opto_[igain].resize(3);
0033   }
0034
0035   for (uint16_t igain = 0; igain < opto_.size(); igain++) {
0036     for (uint16_t ihisto = 0; ihisto < 3; ihisto++) {
0037       // Extra info
0038       std::stringstream extra_info;
0039       extra_info << sistrip::extrainfo::gain_ << igain;
0040       if (ihisto == 0 || ihisto == 1) {
0041         extra_info << sistrip::extrainfo::digital_ << ihisto;
0042       } else {
0043         extra_info << sistrip::extrainfo::baselineRms_;
0044       }
0045
0046       // Title
0047       std::string title = SiStripHistoTitle(sistrip::EXPERT_HISTO,
0048                                             sistrip::OPTO_SCAN,
0049                                             sistrip::FED_KEY,
0050                                             fedKey(),
0051                                             sistrip::LLD_CHAN,
0052                                             connection().lldChannel(),
0053                                             extra_info.str())
0054                               .title();
0055
0056       // Book histo
0057       opto_[igain][ihisto].histo(dqm()->bookProfile(title,
0058                                                     title,
0059                                                     nbins,
0060                                                     0.5,
0061                                                     nbins * 1. + 0.5,  // range is bias setting (1-50)
0062                                                     1024,
0063                                                     -0.5,
0064                                                     1023.5));
0065
0066       opto_[igain][ihisto].vNumOfEntries_.resize(nbins, 0);
0067       opto_[igain][ihisto].vSumOfContents_.resize(nbins, 0);
0068       opto_[igain][ihisto].vSumOfSquares_.resize(nbins, 0);
0069
0070     }  // end loop on histos
0071   }    // end loop on gains
0072 }
0073
0074 // -----------------------------------------------------------------------------
0075 //
0076 void OptoScanTask::fill(const SiStripEventSummary& summary, const edm::DetSet<SiStripRawDigi>& digis) {
0077   //@@ if scope mode length is in trigger fed, then
0078   //@@ can add check here on number of digis
0079   if (digis.data.empty()) {
0080     edm::LogWarning(mlDqmSource_) << "[OptoScanTask::" << __func__ << "]"
0081                                   << " Unexpected number of digis! " << digis.data.size();
0082   } else {
0083     // Retrieve opt bias and gain setting from SiStripEventSummary
0084     uint16_t gain = summary.lldGain();
0085     uint16_t bias = summary.lldBias();
0086
0087     if (gain >= opto_.size()) {
0088       opto_.resize(gain);
0089       for (uint16_t igain = 0; igain < opto_.size(); igain++) {
0090         if (opto_[gain].size() != 3) {
0091           opto_[gain].resize(3);
0092         }
0093       }
0094       edm::LogWarning(mlDqmSource_) << "[OptoScanTask::" << __func__ << "]"
0095                                     << " Unexpected gain value! " << gain;
0096     }
0097
0098     if (bias > 50) {
0099       return;
0100     }  // only use bias settings 1-50
0101
0102     // Find digital "0" and digital "1" levels from tick marks within scope mode data
0103     std::vector<float> baseline;
0104     std::pair<float, float> digital_range;
0105     digital_range.first = sistrip::invalid_;
0106     digital_range.second = sistrip::invalid_;
0107     float baseline_rms = sistrip::invalid_;
0108
0109     locateTicks(digis, digital_range, baseline, baseline_rms);
0110
0111     uint16_t bin = bias - 1;  // fill "bins" (0-49), not bias (1-50)
0112
0113     // Digital "0"
0114     if (digital_range.first < 1. * sistrip::valid_) {
0115       updateHistoSet(opto_[gain][0], bin, digital_range.first);
0116     }
0117
0118     // Digital "1"
0119     if (digital_range.second < 1. * sistrip::valid_) {
0120       updateHistoSet(opto_[gain][1], bin, digital_range.second);
0121     }
0122
0123     // Baseline rms
0124     if (baseline_rms < 1. * sistrip::valid_) {
0125       updateHistoSet(opto_[gain][2], bin, baseline_rms);
0126     }
0127   }
0128 }
0129
0130 // -----------------------------------------------------------------------------
0131 //
0132 void OptoScanTask::update() {
0133   for (uint16_t igain = 0; igain < opto_.size(); igain++) {
0134     for (uint16_t ihisto = 0; ihisto < opto_[igain].size(); ihisto++) {
0135       updateHistoSet(opto_[igain][ihisto]);
0136     }
0137   }
0138 }
0139
0140 // -----------------------------------------------------------------------------
0141 //
0142 void OptoScanTask::locateTicks(const edm::DetSet<SiStripRawDigi>& digis,
0143                                std::pair<float, float>& range,
0144                                std::vector<float>& baseline,
0145                                float& baseline_rms) {
0146   // Copy ADC values and sort
0147   std::vector<uint16_t> adc;
0148   adc.reserve(digis.data.size());
0149   for (uint16_t iadc = 0; iadc < digis.data.size(); iadc++) {
0150     // only take the adc from the first APV (multiplexing alternates the digis)
0151     // this was asked by Karl et al
0152     if (iadc % 2 == 0) {
0153       adc.push_back(digis.data[iadc].adc());
0154     }
0155   }
0156   sort(adc.begin(), adc.end());
0157
0158   // To make sure we have a tickmark, which comes every 70 bxs,
0159   // fully contained in the scope 'frame' we are analyzing
0160   // standard length is 280, sufficient to contain a full apv frame
0161   // in this run there is no frame though, just baseline and tickmarks
0162   if (adc.size() > 70) {
0163     // Define tick mark top" level as "max" ADC values
0164     range.second = adc.back();
0165
0166     // Construct vector to hold "baseline samples" (exclude tick mark samples)
0167     std::vector<uint16_t> truncated;
0168     std::vector<uint16_t>::const_iterator ii = adc.begin();
0169     // remove twice the expected number of tick samples, otherwise you bias the baseline mean and rms
0170     std::vector<uint16_t>::const_iterator jj = adc.end() - 4 * ((adc.size() / 70) + 1);
0171     truncated.resize(jj - ii);
0172     std::copy(ii, jj, truncated.begin());
0173     if (truncated.empty()) {
0174       return;
0175     }
0176
0177     // Calc mean baseline level
0178     float b_mean = 0.;
0179     std::vector<uint16_t>::const_iterator iii = truncated.begin();
0180     std::vector<uint16_t>::const_iterator jjj = truncated.end();
0181     for (; iii != jjj; ++iii) {
0182       b_mean += *iii;
0183     }
0184     b_mean /= (1. * truncated.size());
0185     range.first = b_mean;
0186
0187     // Calc baseline noise
0188     float b_rms = 0.;
0189     std::vector<uint16_t>::const_iterator iiii = truncated.begin();
0190     std::vector<uint16_t>::const_iterator jjjj = truncated.end();
0191     for (; iiii != jjjj; ++iiii) {
0192       b_rms += fabs(*iiii - b_mean);
0193     }
0194     // Set baseline "noise" (requires any possible APV frames are filtered from the data!)
0195     baseline_rms = sqrt(b_rms / (1. * truncated.size()));
0196
0197   } else {
0198     edm::LogWarning(mlDqmSource_) << "[OptoScanTask::" << __func__ << "]"
0199                                   << " Insufficient ADC values: " << adc.size();
0200   }
0201 }