OptoScanAlgorithm.cc

Go to the documentation of this file.

#include "DQM/SiStripCommissioningAnalysis/interface/OptoScanAlgorithm.h"
#include "CondFormats/SiStripObjects/interface/OptoScanAnalysis.h"
#include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
#include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
#include "DQM/SiStripCommissioningAnalysis/src/Utility.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "TProfile.h"
#include "TH1.h"
#include <iostream>
#include <sstream>
#include <iomanip>

using namespace sistrip;

// ----------------------------------------------------------------------------
//
OptoScanAlgorithm::OptoScanAlgorithm(const edm::ParameterSet& pset, OptoScanAnalysis* const anal)
    : CommissioningAlgorithm(anal),
      histos_(4, std::vector<Histo>(3, Histo(nullptr, ""))),
      targetGain_(pset.getParameter<double>("TargetGain")) {
  edm::LogInfo(mlCommissioning_) << "[PedestalsAlgorithm::" << __func__ << "]"
                                 << " Set target gain to: " << targetGain_;
}

// ----------------------------------------------------------------------------
//
void OptoScanAlgorithm::extract(const std::vector<TH1*>& histos) {
  if (!anal()) {
    edm::LogWarning(mlCommissioning_) << "[OptoScanAlgorithm::" << __func__ << "]"
                                      << " NULL pointer to Analysis object!";
    return;
  }

  // Check number of histograms
  if (histos.size() != 12) {
    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::OPTO_SCAN) {
      anal()->addErrorCode(sistrip::unexpectedTask_);
      continue;
    }

    // Extract gain setting and digital high/low info
    uint16_t gain = sistrip::invalid_;
    if (title.extraInfo().find(sistrip::extrainfo::gain_) != std::string::npos) {
      std::stringstream ss;
      ss << title.extraInfo().substr(
          title.extraInfo().find(sistrip::extrainfo::gain_) + (sizeof(sistrip::extrainfo::gain_) - 1), 1);
      ss >> std::dec >> gain;
    }
    uint16_t digital = sistrip::invalid_;
    if (title.extraInfo().find(sistrip::extrainfo::digital_) != std::string::npos) {
      std::stringstream ss;
      ss << title.extraInfo().substr(
          title.extraInfo().find(sistrip::extrainfo::digital_) + (sizeof(sistrip::extrainfo::digital_) - 1), 1);
      ss >> std::dec >> digital;
    }
    bool baseline_rms = false;
    if (title.extraInfo().find(sistrip::extrainfo::baselineRms_) != std::string::npos) {
      baseline_rms = true;
    }

    if (gain <= 3) {
      if (digital <= 1) {
        histos_[gain][digital].first = *ihis;
        histos_[gain][digital].second = (*ihis)->GetName();
      } else if (baseline_rms) {
        histos_[gain][2].first = *ihis;
        histos_[gain][2].second = (*ihis)->GetName();
      } else {
        anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
      }
    } else {
      anal()->addErrorCode(sistrip::unexpectedExtraInfo_);
    }
  }
}

// ----------------------------------------------------------------------------
//
void OptoScanAlgorithm::analyse() {
  if (!anal()) {
    edm::LogWarning(mlCommissioning_) << "[OptoScanAlgorithm::" << __func__ << "]"
                                      << " NULL pointer to base Analysis object!";
    return;
  }

  CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
  OptoScanAnalysis* anal = dynamic_cast<OptoScanAnalysis*>(tmp);
  if (!anal) {
    edm::LogWarning(mlCommissioning_) << "[OptoScanAlgorithm::" << __func__ << "]"
                                      << " NULL pointer to derived Analysis object!";
    return;
  }

  // Iterate through four gain settings
  for (uint16_t igain = 0; igain < 4; igain++) {
    // Select histos appropriate for gain setting
    TH1* base_his = histos_[igain][0].first;
    TH1* peak_his = histos_[igain][1].first;
    TH1* noise_his = histos_[igain][2].first;

    if (!base_his) {
      anal->addErrorCode(sistrip::nullPtr_);
      return;
    }

    if (!peak_his) {
      anal->addErrorCode(sistrip::nullPtr_);
      return;
    }

    if (!noise_his) {
      anal->addErrorCode(sistrip::nullPtr_);
      return;
    }

    TProfile* base_histo = dynamic_cast<TProfile*>(base_his);
    if (!base_histo) {
      anal->addErrorCode(sistrip::nullPtr_);
      return;
    }

    TProfile* peak_histo = dynamic_cast<TProfile*>(peak_his);
    if (!peak_histo) {
      anal->addErrorCode(sistrip::nullPtr_);
      return;
    }

    TProfile* noise_histo = dynamic_cast<TProfile*>(noise_his);
    if (!noise_histo) {
      anal->addErrorCode(sistrip::nullPtr_);
      return;
    }

    // Check histogram binning
    uint16_t nbins = static_cast<uint16_t>(peak_histo->GetNbinsX());
    if (static_cast<uint16_t>(base_histo->GetNbinsX()) != nbins) {
      anal->addErrorCode(sistrip::numberOfBins_);
      if (static_cast<uint16_t>(base_histo->GetNbinsX()) < nbins) {
        nbins = static_cast<uint16_t>(base_histo->GetNbinsX());
      }
    }

    // Some containers
    std::vector<float> peak_contents(0);
    std::vector<float> peak_errors(0);
    std::vector<float> peak_entries(0);
    std::vector<float> base_contents(0);
    std::vector<float> base_errors(0);
    std::vector<float> base_entries(0);
    std::vector<float> noise_contents(0);
    std::vector<float> noise_errors(0);
    std::vector<float> noise_entries(0);
    float peak_max = -1. * sistrip::invalid_;
    float peak_min = 1. * sistrip::invalid_;
    float base_max = -1. * sistrip::invalid_;
    float base_min = 1. * sistrip::invalid_;
    float noise_max = -1. * sistrip::invalid_;
    float noise_min = 1. * sistrip::invalid_;

    // Transfer histogram contents/errors/stats to containers
    for (uint16_t ibin = 0; ibin < nbins; ibin++) {
      // Peak histogram
      peak_contents.push_back(peak_histo->GetBinContent(ibin + 1));
      peak_errors.push_back(peak_histo->GetBinError(ibin + 1));
      peak_entries.push_back(peak_histo->GetBinEntries(ibin + 1));
      if (peak_entries[ibin]) {
        if (peak_contents[ibin] > peak_max) {
          peak_max = peak_contents[ibin];
        }
        if (peak_contents[ibin] < peak_min && ibin) {
          peak_min = peak_contents[ibin];
        }
      }

      // Base histogram
      base_contents.push_back(base_histo->GetBinContent(ibin + 1));
      base_errors.push_back(base_histo->GetBinError(ibin + 1));
      base_entries.push_back(base_histo->GetBinEntries(ibin + 1));
      if (base_entries[ibin]) {
        if (base_contents[ibin] > base_max) {
          base_max = base_contents[ibin];
        }
        if (base_contents[ibin] < base_min && ibin) {
          base_min = base_contents[ibin];
        }
      }

      // Noise histogram
      noise_contents.push_back(noise_histo->GetBinContent(ibin + 1));
      noise_errors.push_back(noise_histo->GetBinError(ibin + 1));
      noise_entries.push_back(noise_histo->GetBinEntries(ibin + 1));
      if (noise_entries[ibin]) {
        if (noise_contents[ibin] > noise_max) {
          noise_max = noise_contents[ibin];
        }
        if (noise_contents[ibin] < noise_min && ibin) {
          noise_min = noise_contents[ibin];
        }
      }
    }

    // Find "zero light" level and error
    //@@ record bias setting used for zero light level
    //@@ zero light error changes wrt gain setting ???
    float zero_light_level = sistrip::invalid_;
    float zero_light_error = sistrip::invalid_;
    for (uint16_t ibin = 0; ibin < nbins; ibin++) {
      if (base_entries[ibin]) {
        zero_light_level = base_contents[ibin];
        zero_light_error = base_errors[ibin];
        break;
      }
    }

    float zero_light_thres = sistrip::invalid_;
    if (zero_light_level <= sistrip::maximum_ && zero_light_error <= sistrip::maximum_) {
      zero_light_thres = zero_light_level + 5. * zero_light_error;
    } else {
      std::stringstream ss;
      ss << sistrip::invalidZeroLightLevel_ << "ForGain" << igain;
      anal->addErrorCode(ss.str());
      continue;
    }

    // Find range of base histogram
    float base_range = base_max - base_min;

    // Find overlapping max/min region that constrains range of linear fit
    float max = peak_max < base_max ? peak_max : base_max;
    float min = peak_min > base_min ? peak_min : base_min;
    float range = max - min;

    // Container identifying whether samples from 'base' histo are above "zero light"
    std::vector<bool> above_zero_light;
    above_zero_light.resize(3, true);

    // Linear fits to top of peak and base curves and one to bottom of base curve
    sistrip::LinearFit peak_high;
    sistrip::LinearFit base_high;
    sistrip::LinearFit base_low;

    // Iterate through histogram bins
    uint16_t peak_bin = 0;
    uint16_t base_bin = 0;
    uint16_t low_bin = 0;
    for (uint16_t ibin = 0; ibin < nbins; ibin++) {
      // Record whether consecutive samples from 'base' histo are above the "zero light" level
      if (base_entries[ibin]) {
        above_zero_light.erase(above_zero_light.begin());
        if (base_contents[ibin] > zero_light_thres) {
          above_zero_light.push_back(true);
        } else {
          above_zero_light.push_back(false);
        }
        if (above_zero_light.size() != 3) {
          above_zero_light.resize(3, false);
        }
      }

      // Linear fit to peak histogram
      if (peak_entries[ibin] && peak_contents[ibin] > (min + 0.2 * range) &&
          peak_contents[ibin] < (min + 0.8 * range)) {
        if (!peak_bin) {
          peak_bin = ibin;
        }
        if ((ibin - peak_bin) < 10) {
          peak_high.add(ibin, peak_contents[ibin]);  //@@ should weight using bin error or bin contents (sqrt(N)/N)
        }
      }
      // Linear fit to base histogram
      if (base_entries[ibin] && base_contents[ibin] > (min + 0.2 * range) &&
          base_contents[ibin] < (min + 0.8 * range)) {
        if (!base_bin) {
          base_bin = ibin;
        }
        if ((ibin - base_bin) < 10) {
          base_high.add(ibin, base_contents[ibin]);  //@@ should weight using bin error or bin contents (sqrt(N)/N)
        }
      }
      // Low linear fit to base histogram
      if (base_entries[ibin] &&
          //@@ above_zero_light[0] && above_zero_light[1] && above_zero_light[2] &&
          base_contents[ibin] > (base_min + 0.2 * base_range) && base_contents[ibin] < (base_min + 0.6 * base_range)) {
        if (!low_bin) {
          low_bin = ibin;
        }
        if ((ibin - low_bin) < 10) {
          base_low.add(ibin, base_contents[ibin]);  //@@ should weight using bin error or bin contents (sqrt(N)/N)
        }
      }
    }

    // Extract width between two curves at midpoint within range
    float mid = min + 0.5 * range;
    sistrip::LinearFit::Params peak_params;
    sistrip::LinearFit::Params base_params;
    peak_high.fit(peak_params);
    base_high.fit(base_params);

    float peak_pos = sistrip::invalid_;
    float base_pos = sistrip::invalid_;
    float width = sistrip::invalid_;
    if (peak_params.b_ > 0.) {
      peak_pos = (mid - peak_params.a_) / peak_params.b_;
    }
    if (base_params.b_ > 0.) {
      base_pos = (mid - base_params.a_) / base_params.b_;
    }
    if (base_pos < sistrip::valid_ && peak_pos < sistrip::valid_) {
      width = base_pos - peak_pos;
    }

    // Extrapolate to zero light to find "lift off"
    sistrip::LinearFit::Params low_params;
    base_low.fit(low_params);
    float lift_off = sistrip::invalid_;
    if (low_params.b_ > 0.) {
      lift_off = (zero_light_level - low_params.a_) / low_params.b_;
    }

    // ---------- Set all parameters ----------

    // Slope of baseline
    if (low_params.b_ > 0.) {
      anal->baseSlope_[igain] = low_params.b_;
    }

    // Check "lift off" value and set bias setting accordingly
    if (lift_off <= sistrip::maximum_) {
      anal->bias_[igain] = static_cast<uint16_t>(lift_off) + 2;
    } else {
      anal->bias_[igain] = OptoScanAnalysis::defaultBiasSetting_;
    }

    // Calculate "lift off" (in mA)
    if (lift_off <= sistrip::maximum_) {
      anal->liftOff_[igain] = 0.45 * lift_off;
    }

    // Calculate laser threshold (in mA)
    if (width < sistrip::invalid_) {
      anal->threshold_[igain] = 0.45 * (lift_off - width / 2.);
    }

    // Set "zero light" level
    anal->zeroLight_[igain] = zero_light_level;

    // Set link noise
    uint16_t bin_number = sistrip::invalid_;
    if (anal->threshold_[igain] < sistrip::valid_) {
      // Old calculation, used in commissioning in 2008
      //   always leads to zero link noise
      //   bin_number = static_cast<uint16_t>( anal->threshold_[igain] / 0.45 );
      // New calculation asked by Karl et al, for commissioning in 2009
      bin_number = (uint16_t)(lift_off + width / 3.);
    }
    if (bin_number < sistrip::valid_) {
      if (bin_number < noise_contents.size()) {
        anal->linkNoise_[igain] = noise_contents[bin_number];
      } else {
        anal->addErrorCode(sistrip::unexpectedBinNumber_);
      }
    }

    // Calculate tick mark height
    if (low_params.b_ <= sistrip::maximum_ && width <= sistrip::maximum_) {
      anal->tickHeight_[igain] = width * low_params.b_;
    }

    // Set measured gain
    if (anal->tickHeight_[igain] < sistrip::invalid_ - 1.) {
      anal->measGain_[igain] = anal->tickHeight_[igain] * OptoScanAnalysis::fedAdcGain_ / 0.800;
    } else {
      anal->measGain_[igain] = sistrip::invalid_;
    }

  }  // gain loop

  // Iterate through four gain settings and identify optimum gain setting
  const float target_gain =
      targetGain_;  //0.863; // changed from 0.8 to avoid choice of low tickheights (Xtof, SL, 15/6/2009)

  float diff_in_gain = sistrip::invalid_;
  for (uint16_t igain = 0; igain < 4; igain++) {
    // Check for sensible gain value
    if (anal->measGain_[igain] > sistrip::maximum_) {
      continue;
    }

    // Find optimum gain setting
    if (fabs(anal->measGain_[igain] - target_gain) < diff_in_gain) {
      anal->gain_ = igain;
      diff_in_gain = fabs(anal->measGain_[igain] - target_gain);
    }
  }

  // Check optimum gain setting
  if (anal->gain_ > sistrip::maximum_) {
    anal->gain_ = OptoScanAnalysis::defaultGainSetting_;
  }
}

// ----------------------------------------------------------------------------
//
CommissioningAlgorithm::Histo OptoScanAlgorithm::histo(const uint16_t& gain, const uint16_t& digital_level) const {
  if (gain <= 3 && digital_level <= 1) {
    return histos_[gain][digital_level];
  } else {
    return Histo(nullptr, "");
  }
}