SiStripCommissioningAnalysis/src/OptoScanAlgorithm.cc

0001 #include "DQM/SiStripCommissioningAnalysis/interface/OptoScanAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/OptoScanAnalysis.h"
0003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0004 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
0005 #include "DQM/SiStripCommissioningAnalysis/src/Utility.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "TProfile.h"
0008 #include "TH1.h"
0009 #include <iostream>
0010 #include <sstream>
0011 #include <iomanip>
0012
0013 using namespace sistrip;
0014
0015 // ----------------------------------------------------------------------------
0016 //
0017 OptoScanAlgorithm::OptoScanAlgorithm(const edm::ParameterSet& pset, OptoScanAnalysis* const anal)
0018     : CommissioningAlgorithm(anal),
0019       histos_(4, std::vector<Histo>(3, Histo(nullptr, ""))),
0020       targetGain_(pset.getParameter<double>("TargetGain")) {
0021   edm::LogInfo(mlCommissioning_) << "[PedestalsAlgorithm::" << __func__ << "]"
0022                                  << " Set target gain to: " << targetGain_;
0023 }
0024
0025 // ----------------------------------------------------------------------------
0026 //
0027 void OptoScanAlgorithm::extract(const std::vector<TH1*>& histos) {
0028   if (!anal()) {
0029     edm::LogWarning(mlCommissioning_) << "[OptoScanAlgorithm::" << __func__ << "]"
0030                                       << " NULL pointer to Analysis object!";
0031     return;
0032   }
0033
0034   // Check number of histograms
0035   if (histos.size() != 12) {
0036     anal()->addErrorCode(sistrip::numberOfHistos_);
0037   }
0038
0039   // Extract FED key from histo title
0040   if (!histos.empty())
0041     anal()->fedKey(extractFedKey(histos.front()));
0042
0043   // Extract histograms
0044   std::vector<TH1*>::const_iterator ihis = histos.begin();
0045   for (; ihis != histos.end(); ihis++) {
0046     // Check for NULL pointer
0047     if (!(*ihis)) {
0048       continue;
0049     }
0050
0051     // Check name
0052     SiStripHistoTitle title((*ihis)->GetName());
0053     if (title.runType() != sistrip::OPTO_SCAN) {
0054       anal()->addErrorCode(sistrip::unexpectedTask_);
0055       continue;
0056     }
0057
0058     // Extract gain setting and digital high/low info
0059     uint16_t gain = sistrip::invalid_;
0060     if (title.extraInfo().find(sistrip::extrainfo::gain_) != std::string::npos) {
0061       std::stringstream ss;
0062       ss << title.extraInfo().substr(
0063           title.extraInfo().find(sistrip::extrainfo::gain_) + (sizeof(sistrip::extrainfo::gain_) - 1), 1);
0064       ss >> std::dec >> gain;
0065     }
0066     uint16_t digital = sistrip::invalid_;
0067     if (title.extraInfo().find(sistrip::extrainfo::digital_) != std::string::npos) {
0068       std::stringstream ss;
0069       ss << title.extraInfo().substr(
0070           title.extraInfo().find(sistrip::extrainfo::digital_) + (sizeof(sistrip::extrainfo::digital_) - 1), 1);
0071       ss >> std::dec >> digital;
0072     }
0073     bool baseline_rms = false;
0074     if (title.extraInfo().find(sistrip::extrainfo::baselineRms_) != std::string::npos) {
0075       baseline_rms = true;
0076     }
0077
0078     if (gain <= 3) {
0079       if (digital <= 1) {
0080         histos_[gain][digital].first = *ihis;
0081         histos_[gain][digital].second = (*ihis)->GetName();
0082       } else if (baseline_rms) {
0083         histos_[gain][2].first = *ihis;
0084         histos_[gain][2].second = (*ihis)->GetName();
0085       } else {
0086         anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
0087       }
0088     } else {
0089       anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
0090     }
0091   }
0092 }
0093
0094 // ----------------------------------------------------------------------------
0095 //
0096 void OptoScanAlgorithm::analyse() {
0097   if (!anal()) {
0098     edm::LogWarning(mlCommissioning_) << "[OptoScanAlgorithm::" << __func__ << "]"
0099                                       << " NULL pointer to base Analysis object!";
0100     return;
0101   }
0102
0103   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0104   OptoScanAnalysis* anal = dynamic_cast<OptoScanAnalysis*>(tmp);
0105   if (!anal) {
0106     edm::LogWarning(mlCommissioning_) << "[OptoScanAlgorithm::" << __func__ << "]"
0107                                       << " NULL pointer to derived Analysis object!";
0108     return;
0109   }
0110
0111   // Iterate through four gain settings
0112   for (uint16_t igain = 0; igain < 4; igain++) {
0113     // Select histos appropriate for gain setting
0114     TH1* base_his = histos_[igain][0].first;
0115     TH1* peak_his = histos_[igain][1].first;
0116     TH1* noise_his = histos_[igain][2].first;
0117
0118     if (!base_his) {
0119       anal->addErrorCode(sistrip::nullPtr_);
0120       return;
0121     }
0122
0123     if (!peak_his) {
0124       anal->addErrorCode(sistrip::nullPtr_);
0125       return;
0126     }
0127
0128     if (!noise_his) {
0129       anal->addErrorCode(sistrip::nullPtr_);
0130       return;
0131     }
0132
0133     TProfile* base_histo = dynamic_cast<TProfile*>(base_his);
0134     if (!base_histo) {
0135       anal->addErrorCode(sistrip::nullPtr_);
0136       return;
0137     }
0138
0139     TProfile* peak_histo = dynamic_cast<TProfile*>(peak_his);
0140     if (!peak_histo) {
0141       anal->addErrorCode(sistrip::nullPtr_);
0142       return;
0143     }
0144
0145     TProfile* noise_histo = dynamic_cast<TProfile*>(noise_his);
0146     if (!noise_histo) {
0147       anal->addErrorCode(sistrip::nullPtr_);
0148       return;
0149     }
0150
0151     // Check histogram binning
0152     uint16_t nbins = static_cast<uint16_t>(peak_histo->GetNbinsX());
0153     if (static_cast<uint16_t>(base_histo->GetNbinsX()) != nbins) {
0154       anal->addErrorCode(sistrip::numberOfBins_);
0155       if (static_cast<uint16_t>(base_histo->GetNbinsX()) < nbins) {
0156         nbins = static_cast<uint16_t>(base_histo->GetNbinsX());
0157       }
0158     }
0159
0160     // Some containers
0161     std::vector<float> peak_contents(0);
0162     std::vector<float> peak_errors(0);
0163     std::vector<float> peak_entries(0);
0164     std::vector<float> base_contents(0);
0165     std::vector<float> base_errors(0);
0166     std::vector<float> base_entries(0);
0167     std::vector<float> noise_contents(0);
0168     std::vector<float> noise_errors(0);
0169     std::vector<float> noise_entries(0);
0170     float peak_max = -1. * sistrip::invalid_;
0171     float peak_min = 1. * sistrip::invalid_;
0172     float base_max = -1. * sistrip::invalid_;
0173     float base_min = 1. * sistrip::invalid_;
0174     float noise_max = -1. * sistrip::invalid_;
0175     float noise_min = 1. * sistrip::invalid_;
0176
0177     // Transfer histogram contents/errors/stats to containers
0178     for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0179       // Peak histogram
0180       peak_contents.push_back(peak_histo->GetBinContent(ibin + 1));
0181       peak_errors.push_back(peak_histo->GetBinError(ibin + 1));
0182       peak_entries.push_back(peak_histo->GetBinEntries(ibin + 1));
0183       if (peak_entries[ibin]) {
0184         if (peak_contents[ibin] > peak_max) {
0185           peak_max = peak_contents[ibin];
0186         }
0187         if (peak_contents[ibin] < peak_min && ibin) {
0188           peak_min = peak_contents[ibin];
0189         }
0190       }
0191
0192       // Base histogram
0193       base_contents.push_back(base_histo->GetBinContent(ibin + 1));
0194       base_errors.push_back(base_histo->GetBinError(ibin + 1));
0195       base_entries.push_back(base_histo->GetBinEntries(ibin + 1));
0196       if (base_entries[ibin]) {
0197         if (base_contents[ibin] > base_max) {
0198           base_max = base_contents[ibin];
0199         }
0200         if (base_contents[ibin] < base_min && ibin) {
0201           base_min = base_contents[ibin];
0202         }
0203       }
0204
0205       // Noise histogram
0206       noise_contents.push_back(noise_histo->GetBinContent(ibin + 1));
0207       noise_errors.push_back(noise_histo->GetBinError(ibin + 1));
0208       noise_entries.push_back(noise_histo->GetBinEntries(ibin + 1));
0209       if (noise_entries[ibin]) {
0210         if (noise_contents[ibin] > noise_max) {
0211           noise_max = noise_contents[ibin];
0212         }
0213         if (noise_contents[ibin] < noise_min && ibin) {
0214           noise_min = noise_contents[ibin];
0215         }
0216       }
0217     }
0218
0219     // Find "zero light" level and error
0220     //@@ record bias setting used for zero light level
0221     //@@ zero light error changes wrt gain setting ???
0222     float zero_light_level = sistrip::invalid_;
0223     float zero_light_error = sistrip::invalid_;
0224     for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0225       if (base_entries[ibin]) {
0226         zero_light_level = base_contents[ibin];
0227         zero_light_error = base_errors[ibin];
0228         break;
0229       }
0230     }
0231
0232     float zero_light_thres = sistrip::invalid_;
0233     if (zero_light_level <= sistrip::maximum_ && zero_light_error <= sistrip::maximum_) {
0234       zero_light_thres = zero_light_level + 5. * zero_light_error;
0235     } else {
0236       std::stringstream ss;
0237       ss << sistrip::invalidZeroLightLevel_ << "ForGain" << igain;
0238       anal->addErrorCode(ss.str());
0239       continue;
0240     }
0241
0242     // Find range of base histogram
0243     float base_range = base_max - base_min;
0244
0245     // Find overlapping max/min region that constrains range of linear fit
0246     float max = peak_max < base_max ? peak_max : base_max;
0247     float min = peak_min > base_min ? peak_min : base_min;
0248     float range = max - min;
0249
0250     // Container identifying whether samples from 'base' histo are above "zero light"
0251     std::vector<bool> above_zero_light;
0252     above_zero_light.resize(3, true);
0253
0254     // Linear fits to top of peak and base curves and one to bottom of base curve
0255     sistrip::LinearFit peak_high;
0256     sistrip::LinearFit base_high;
0257     sistrip::LinearFit base_low;
0258
0259     // Iterate through histogram bins
0260     uint16_t peak_bin = 0;
0261     uint16_t base_bin = 0;
0262     uint16_t low_bin = 0;
0263     for (uint16_t ibin = 0; ibin < nbins; ibin++) {
0264       // Record whether consecutive samples from 'base' histo are above the "zero light" level
0265       if (base_entries[ibin]) {
0266         above_zero_light.erase(above_zero_light.begin());
0267         if (base_contents[ibin] > zero_light_thres) {
0268           above_zero_light.push_back(true);
0269         } else {
0270           above_zero_light.push_back(false);
0271         }
0272         if (above_zero_light.size() != 3) {
0273           above_zero_light.resize(3, false);
0274         }
0275       }
0276
0277       // Linear fit to peak histogram
0278       if (peak_entries[ibin] && peak_contents[ibin] > (min + 0.2 * range) &&
0279           peak_contents[ibin] < (min + 0.8 * range)) {
0280         if (!peak_bin) {
0281           peak_bin = ibin;
0282         }
0283         if ((ibin - peak_bin) < 10) {
0284           peak_high.add(ibin, peak_contents[ibin]);  //@@ should weight using bin error or bin contents (sqrt(N)/N)
0285         }
0286       }
0287       // Linear fit to base histogram
0288       if (base_entries[ibin] && base_contents[ibin] > (min + 0.2 * range) &&
0289           base_contents[ibin] < (min + 0.8 * range)) {
0290         if (!base_bin) {
0291           base_bin = ibin;
0292         }
0293         if ((ibin - base_bin) < 10) {
0294           base_high.add(ibin, base_contents[ibin]);  //@@ should weight using bin error or bin contents (sqrt(N)/N)
0295         }
0296       }
0297       // Low linear fit to base histogram
0298       if (base_entries[ibin] &&
0299           //@@ above_zero_light[0] && above_zero_light[1] && above_zero_light[2] &&
0300           base_contents[ibin] > (base_min + 0.2 * base_range) && base_contents[ibin] < (base_min + 0.6 * base_range)) {
0301         if (!low_bin) {
0302           low_bin = ibin;
0303         }
0304         if ((ibin - low_bin) < 10) {
0305           base_low.add(ibin, base_contents[ibin]);  //@@ should weight using bin error or bin contents (sqrt(N)/N)
0306         }
0307       }
0308     }
0309
0310     // Extract width between two curves at midpoint within range
0311     float mid = min + 0.5 * range;
0312     sistrip::LinearFit::Params peak_params;
0313     sistrip::LinearFit::Params base_params;
0314     peak_high.fit(peak_params);
0315     base_high.fit(base_params);
0316
0317     float peak_pos = sistrip::invalid_;
0318     float base_pos = sistrip::invalid_;
0319     float width = sistrip::invalid_;
0320     if (peak_params.b_ > 0.) {
0321       peak_pos = (mid - peak_params.a_) / peak_params.b_;
0322     }
0323     if (base_params.b_ > 0.) {
0324       base_pos = (mid - base_params.a_) / base_params.b_;
0325     }
0326     if (base_pos < sistrip::valid_ && peak_pos < sistrip::valid_) {
0327       width = base_pos - peak_pos;
0328     }
0329
0330     // Extrapolate to zero light to find "lift off"
0331     sistrip::LinearFit::Params low_params;
0332     base_low.fit(low_params);
0333     float lift_off = sistrip::invalid_;
0334     if (low_params.b_ > 0.) {
0335       lift_off = (zero_light_level - low_params.a_) / low_params.b_;
0336     }
0337
0338     // ---------- Set all parameters ----------
0339
0340     // Slope of baseline
0341     if (low_params.b_ > 0.) {
0342       anal->baseSlope_[igain] = low_params.b_;
0343     }
0344
0345     // Check "lift off" value and set bias setting accordingly
0346     if (lift_off <= sistrip::maximum_) {
0347       anal->bias_[igain] = static_cast<uint16_t>(lift_off) + 2;
0348     } else {
0349       anal->bias_[igain] = OptoScanAnalysis::defaultBiasSetting_;
0350     }
0351
0352     // Calculate "lift off" (in mA)
0353     if (lift_off <= sistrip::maximum_) {
0354       anal->liftOff_[igain] = 0.45 * lift_off;
0355     }
0356
0357     // Calculate laser threshold (in mA)
0358     if (width < sistrip::invalid_) {
0359       anal->threshold_[igain] = 0.45 * (lift_off - width / 2.);
0360     }
0361
0362     // Set "zero light" level
0363     anal->zeroLight_[igain] = zero_light_level;
0364
0365     // Set link noise
0366     uint16_t bin_number = sistrip::invalid_;
0367     if (anal->threshold_[igain] < sistrip::valid_) {
0368       // Old calculation, used in commissioning in 2008
0369       //   always leads to zero link noise
0370       //   bin_number = static_cast<uint16_t>( anal->threshold_[igain] / 0.45 );
0371       // New calculation asked by Karl et al, for commissioning in 2009
0372       bin_number = (uint16_t)(lift_off + width / 3.);
0373     }
0374     if (bin_number < sistrip::valid_) {
0375       if (bin_number < noise_contents.size()) {
0376         anal->linkNoise_[igain] = noise_contents[bin_number];
0377       } else {
0378         anal->addErrorCode(sistrip::unexpectedBinNumber_);
0379       }
0380     }
0381
0382     // Calculate tick mark height
0383     if (low_params.b_ <= sistrip::maximum_ && width <= sistrip::maximum_) {
0384       anal->tickHeight_[igain] = width * low_params.b_;
0385     }
0386
0387     // Set measured gain
0388     if (anal->tickHeight_[igain] < sistrip::invalid_ - 1.) {
0389       anal->measGain_[igain] = anal->tickHeight_[igain] * OptoScanAnalysis::fedAdcGain_ / 0.800;
0390     } else {
0391       anal->measGain_[igain] = sistrip::invalid_;
0392     }
0393
0394   }  // gain loop
0395
0396   // Iterate through four gain settings and identify optimum gain setting
0397   const float target_gain =
0398       targetGain_;  //0.863; // changed from 0.8 to avoid choice of low tickheights (Xtof, SL, 15/6/2009)
0399
0400   float diff_in_gain = sistrip::invalid_;
0401   for (uint16_t igain = 0; igain < 4; igain++) {
0402     // Check for sensible gain value
0403     if (anal->measGain_[igain] > sistrip::maximum_) {
0404       continue;
0405     }
0406
0407     // Find optimum gain setting
0408     if (fabs(anal->measGain_[igain] - target_gain) < diff_in_gain) {
0409       anal->gain_ = igain;
0410       diff_in_gain = fabs(anal->measGain_[igain] - target_gain);
0411     }
0412   }
0413
0414   // Check optimum gain setting
0415   if (anal->gain_ > sistrip::maximum_) {
0416     anal->gain_ = OptoScanAnalysis::defaultGainSetting_;
0417   }
0418 }
0419
0420 // ----------------------------------------------------------------------------
0421 //
0422 CommissioningAlgorithm::Histo OptoScanAlgorithm::histo(const uint16_t& gain, const uint16_t& digital_level) const {
0423   if (gain <= 3 && digital_level <= 1) {
0424     return histos_[gain][digital_level];
0425   } else {
0426     return Histo(nullptr, "");
0427   }
0428 }