CSCConditions.cc

Go to the documentation of this file.

#include "CalibMuon/CSCCalibration/interface/CSCConditions.h"

#include "CalibMuon/CSCCalibration/interface/CSCChannelMapperBase.h"
#include "CalibMuon/CSCCalibration/interface/CSCIndexerBase.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/ESWatcher.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "CondFormats/CSCObjects/interface/CSCChamberTimeCorrections.h"
#include "CondFormats/CSCObjects/interface/CSCDBChipSpeedCorrection.h"
#include "CondFormats/CSCObjects/interface/CSCDBCrosstalk.h"
#include "CondFormats/CSCObjects/interface/CSCDBGains.h"
#include "CondFormats/CSCObjects/interface/CSCDBGasGainCorrection.h"
#include "CondFormats/CSCObjects/interface/CSCDBPedestals.h"

#include "CondFormats/CSCObjects/interface/CSCBadChambers.h"
#include "CondFormats/CSCObjects/interface/CSCBadStrips.h"
#include "CondFormats/CSCObjects/interface/CSCBadWires.h"

CSCConditions::CSCConditions(const edm::ParameterSet &ps, edm::ConsumesCollector cc)
    : theGains(),
      theCrosstalk(),
      thePedestals(),
      theNoiseMatrix(),
      theBadStrips(),
      theBadWires(),
      theBadChambers(),
      theChipCorrections(),
      theChamberTimingCorrections(),
      theGasGainCorrections(),
      indexer_(nullptr),
      mapper_(nullptr),
      readBadChannels_(false),
      readBadChambers_(false),
      useTimingCorrections_(false),
      useGasGainCorrections_(false),
      idOfBadChannelWords_(CSCDetId()),
      badStripWord_(0),
      badWireWord_(0),
      theAverageGain(-1.0) {
  readBadChannels_ = ps.getParameter<bool>("readBadChannels");
  readBadChambers_ = ps.getParameter<bool>("readBadChambers");
  useTimingCorrections_ = ps.getParameter<bool>("CSCUseTimingCorrections");
  useGasGainCorrections_ = ps.getParameter<bool>("CSCUseGasGainCorrections");
  indexerToken_ = cc.esConsumes<CSCIndexerBase, CSCIndexerRecord>();
  mapperToken_ = cc.esConsumes<CSCChannelMapperBase, CSCChannelMapperRecord>();
  gainsToken_ = cc.esConsumes<CSCDBGains, CSCDBGainsRcd>();
  crosstalkToken_ = cc.esConsumes<CSCDBCrosstalk, CSCDBCrosstalkRcd>();
  pedestalsToken_ = cc.esConsumes<CSCDBPedestals, CSCDBPedestalsRcd>();
  noiseMatrixToken_ = cc.esConsumes<CSCDBNoiseMatrix, CSCDBNoiseMatrixRcd>();
  if (useTimingCorrections()) {
    chipCorrectionsToken_ = cc.esConsumes<CSCDBChipSpeedCorrection, CSCDBChipSpeedCorrectionRcd>();
    chamberTimingCorrectionsToken_ = cc.esConsumes<CSCChamberTimeCorrections, CSCChamberTimeCorrectionsRcd>();
  }
  if (readBadChannels()) {
    badStripsToken_ = cc.esConsumes<CSCBadStrips, CSCBadStripsRcd>();
    badWiresToken_ = cc.esConsumes<CSCBadWires, CSCBadWiresRcd>();
  }
  if (readBadChambers()) {
    badChambersToken_ = cc.esConsumes<CSCBadChambers, CSCBadChambersRcd>();
  }
  if (useGasGainCorrections()) {
    gasGainCorrectionsToken_ = cc.esConsumes<CSCDBGasGainCorrection, CSCDBGasGainCorrectionRcd>();
  }
}

CSCConditions::~CSCConditions() {}

void CSCConditions::initializeEvent(const edm::EventSetup &es) {
  // Algorithms
  indexer_ = es.getHandle(indexerToken_);
  mapper_ = es.getHandle(mapperToken_);

  // Strip gains
  theGains = es.getHandle(gainsToken_);
  // Strip X-talk
  theCrosstalk = es.getHandle(crosstalkToken_);
  // Strip pedestals
  thePedestals = es.getHandle(pedestalsToken_);
  // Strip autocorrelation noise matrix
  theNoiseMatrix = es.getHandle(noiseMatrixToken_);

  if (useTimingCorrections()) {
    // Buckeye chip speeds
    theChipCorrections = es.getHandle(chipCorrectionsToken_);
    // Cable lengths from chambers to peripheral crate and additional chamber
    // level timing correction
    theChamberTimingCorrections = es.getHandle(chamberTimingCorrectionsToken_);
  }

  if (readBadChannels()) {
    // Bad strip channels
    theBadStrips = es.getHandle(badStripsToken_);
    // Bad wiregroup channels
    theBadWires = es.getHandle(badWiresToken_);

    //@@    if( badStripsWatcher_.check( es ) ) {
    //      fillBadStripWords();
    //@@    }
    //@@    if( badWiresWatcher_.check( es ) ) {
    //      fillBadWireWords();
    //@    }
  }

  // Has GainsRcd changed?
  if (gainsWatcher_.check(es)) {  // Yes...
    theAverageGain = -1.0;        // ...reset, so next access will recalculate it
  }

  if (readBadChambers()) {
    // Entire bad chambers
    theBadChambers = es.getHandle(badChambersToken_);
  }

  if (useGasGainCorrections()) {
    theGasGainCorrections = es.getHandle(gasGainCorrectionsToken_);
  }

  //  print();
}

void CSCConditions::fillBadChannelWords(const CSCDetId &id) {
  // input CSCDetId is expected to be an offline value i.e. different for ME1/1A
  // and ME1/1B

  // Only update content if necessary
  if (id != idOfBadChannelWords()) {
    // store offline CSCDetId for the two bad channel words
    setIdOfBadChannelWords(id);

    // reset to all zeroes
    badStripWord_.reset();
    badWireWord_.reset();

    if (readBadChannels()) {
      // convert to online CSCDetId since that is how conditions data are stored
      CSCDetId idraw = mapper_->rawCSCDetId(id);
      fillBadStripWord(idraw);
      fillBadWireWord(idraw);
    }
  }
}


void CSCConditions::fillBadStripWord(const CSCDetId &id) {
  // Input CSCDetId is expected to be a 'raw' value

  // Find linear index of chamber for input CSCDetId
  int inputIndex = indexer_->chamberIndex(id);
  short inputLayer = id.layer();

  // Does this chamber occur in bad channel list? If so, unpack its bad channels

  // chambers is a vector<BadChamber>
  // channels is a vector<BadChannel>
  // Each BadChamber contains its index (1-468 or 540 w. ME42), the no. of bad
  // channels, and the index within vector<BadChannel> where this chamber's bad
  // channels start.

  for (size_t i = 0; i < theBadStrips->chambers.size(); ++i) {  // loop over bad chambers
    int indexc = theBadStrips->chambers[i].chamber_index;
    if (indexc != inputIndex)
      continue;  // next iteration if not a match

    int start = theBadStrips->chambers[i].pointer;
    int nbad = theBadStrips->chambers[i].bad_channels;

    for (int j = start - 1; j < start - 1 + nbad; ++j) {  // bad channels in this chamber
      short lay = theBadStrips->channels[j].layer;        // value 1-6
      if (lay != inputLayer)
        continue;

      short chan = theBadStrips->channels[j].channel;  // value 1-80 (->112 for unganged ME1/1A)
      // Flags so far unused (and unset in conditins data)
      //    short f1 = theBadStrips->channels[j].flag1;
      //    short f2 = theBadStrips->channels[j].flag2;
      //    short f3 = theBadStrips->channels[j].flag3;
      badStripWord_.set(chan - 1, true);  // set bit 0-79 (111) in 80 (112)-bit
                                          // bitset representing this layer
    }                                     // j
  }                                       // i
}

void CSCConditions::fillBadWireWord(const CSCDetId &id) {
  // Input CSCDetId is expected to be a 'raw' value

  // Find linear index of chamber for input CSCDetId
  int inputIndex = indexer_->chamberIndex(id);
  short inputLayer = id.layer();

  // unpack what we've read from theBadWires

  for (size_t i = 0; i < theBadWires->chambers.size(); ++i) {  // loop over bad chambers
    int indexc = theBadWires->chambers[i].chamber_index;

    if (indexc != inputIndex)
      continue;  // next iteration if not a match

    int start = theBadWires->chambers[i].pointer;
    int nbad = theBadWires->chambers[i].bad_channels;

    for (int j = start - 1; j < start - 1 + nbad; ++j) {  // bad channels in this chamber
      short lay = theBadWires->channels[j].layer;         // value 1-6
      if (lay != inputLayer)
        continue;

      short chan = theBadWires->channels[j].channel;  // value 1-112
      //    short f1 = theBadWires->channels[j].flag1;
      //    short f2 = theBadWires->channels[j].flag2;
      //    short f3 = theBadWires->channels[j].flag3;
      badWireWord_.set(chan - 1,
                       true);  // set bit 0-111 in 112-bit bitset representing this layer
    }                          // j
  }                            // i
}

bool CSCConditions::isInBadChamber(const CSCDetId &id) const {
  //@@ We do not consider the possibility of having ME1/1A & ME1/1B
  // independently 'bad'.
  //@@ To do that we would need to define separate chamber indexes for ME1/1A &
  // ME1/1B.

  if (readBadChambers()) {
    CSCDetId idraw = mapper_->rawCSCDetId(id);
    int index = indexer_->chamberIndex(idraw);
    return theBadChambers->isInBadChamber(index);
  } else
    return false;
}

float CSCConditions::gain(const CSCDetId &id, int geomChannel) const {
  assert(theGains.isValid());
  CSCDetId idraw = mapper_->rawCSCDetId(id);
  int iraw = mapper_->rawStripChannel(id, geomChannel);
  int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
  return float(theGains->gain(index)) / theGains->scale();
}

float CSCConditions::pedestal(const CSCDetId &id, int geomChannel) const {
  assert(thePedestals.isValid());
  CSCDetId idraw = mapper_->rawCSCDetId(id);
  int iraw = mapper_->rawStripChannel(id, geomChannel);
  int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
  return float(thePedestals->pedestal(index)) / thePedestals->scale_ped();
}

float CSCConditions::pedestalSigma(const CSCDetId &id, int geomChannel) const {
  assert(thePedestals.isValid());
  CSCDetId idraw = mapper_->rawCSCDetId(id);
  int iraw = mapper_->rawStripChannel(id, geomChannel);
  int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
  return float(thePedestals->pedestal_rms(index)) / thePedestals->scale_rms();
}

float CSCConditions::crosstalkIntercept(const CSCDetId &id, int geomChannel, bool leftRight) const {
  assert(theCrosstalk.isValid());
  CSCDetId idraw = mapper_->rawCSCDetId(id);
  int iraw = mapper_->rawStripChannel(id, geomChannel);
  int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
  // resistive fraction is at the peak, where t=0
  return leftRight ? float(theCrosstalk->rinter(index)) / theCrosstalk->iscale()
                   : float(theCrosstalk->linter(index)) / theCrosstalk->iscale();
}

float CSCConditions::crosstalkSlope(const CSCDetId &id, int geomChannel, bool leftRight) const {
  assert(theCrosstalk.isValid());
  CSCDetId idraw = mapper_->rawCSCDetId(id);
  int iraw = mapper_->rawStripChannel(id, geomChannel);
  int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
  // resistive fraction is at the peak, where t=0
  return leftRight ? float(theCrosstalk->rslope(index)) / theCrosstalk->sscale()
                   : float(theCrosstalk->lslope(index)) / theCrosstalk->sscale();
}

const CSCDBNoiseMatrix::Item &CSCConditions::noiseMatrix(const CSCDetId &id, int geomChannel) const {
  //@@ BEWARE - THIS FUNCTION DOES NOT APPLY SCALE FACTOR USED IN PACKING VALUES
  // IN CONDITIONS DATA
  //@@ MAY BE AN ERROR? WHO WOULD WANT ACCESS WITHOUT IT?

  assert(theNoiseMatrix.isValid());
  CSCDetId idraw = mapper_->rawCSCDetId(id);
  int iraw = mapper_->rawStripChannel(id, geomChannel);
  int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
  return theNoiseMatrix->item(index);
}

void CSCConditions::noiseMatrixElements(const CSCDetId &id, int geomChannel, std::vector<float> &me) const {
  assert(me.size() > 11);
  const CSCDBNoiseMatrix::Item &item = noiseMatrix(id, geomChannel);  // i.e. the function above
  me[0] = float(item.elem33) / theNoiseMatrix->scale();
  me[1] = float(item.elem34) / theNoiseMatrix->scale();
  me[2] = float(item.elem35) / theNoiseMatrix->scale();
  me[3] = float(item.elem44) / theNoiseMatrix->scale();
  me[4] = float(item.elem45) / theNoiseMatrix->scale();
  me[5] = float(item.elem46) / theNoiseMatrix->scale();
  me[6] = float(item.elem55) / theNoiseMatrix->scale();
  me[7] = float(item.elem56) / theNoiseMatrix->scale();
  me[8] = float(item.elem57) / theNoiseMatrix->scale();
  me[9] = float(item.elem66) / theNoiseMatrix->scale();
  me[10] = float(item.elem67) / theNoiseMatrix->scale();
  me[11] = float(item.elem77) / theNoiseMatrix->scale();
}

void CSCConditions::crossTalk(const CSCDetId &id, int geomChannel, std::vector<float> &ct) const {
  assert(theCrosstalk.isValid());
  CSCDetId idraw = mapper_->rawCSCDetId(id);
  int iraw = mapper_->rawStripChannel(id, geomChannel);
  int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!

  ct[0] = float(theCrosstalk->lslope(index)) / theCrosstalk->sscale();
  ct[1] = float(theCrosstalk->linter(index)) / theCrosstalk->iscale();
  ct[2] = float(theCrosstalk->rslope(index)) / theCrosstalk->sscale();
  ct[3] = float(theCrosstalk->rinter(index)) / theCrosstalk->iscale();
}

float CSCConditions::chipCorrection(const CSCDetId &id, int geomChannel) const {
  if (useTimingCorrections()) {
    assert(theChipCorrections.isValid());
    CSCDetId idraw = mapper_->rawCSCDetId(id);
    int iraw = mapper_->rawStripChannel(id, geomChannel);
    int ichip = indexer_->chipIndex(iraw);              // converts 1-80 to 1-5 (chip#, CFEB#)
    int index = indexer_->chipIndex(idraw, ichip) - 1;  // NOTE THE MINUS ONE!
    return float(theChipCorrections->value(index)) / theChipCorrections->scale();
  } else
    return 0;
}
float CSCConditions::chamberTimingCorrection(const CSCDetId &id) const {
  if (useTimingCorrections()) {
    assert(theChamberTimingCorrections.isValid());
    CSCDetId idraw = mapper_->rawCSCDetId(id);
    int index = indexer_->chamberIndex(idraw) - 1;  // NOTE THE MINUS ONE!
    return float(
        theChamberTimingCorrections->item(index).cfeb_tmb_skew_delay * 1. / theChamberTimingCorrections->precision() +
        theChamberTimingCorrections->item(index).cfeb_timing_corr * 1. / theChamberTimingCorrections->precision() +
        (theChamberTimingCorrections->item(index).cfeb_cable_delay * 25.));
  } else
    return 0;
}
float CSCConditions::anodeBXoffset(const CSCDetId &id) const {
  if (useTimingCorrections()) {
    assert(theChamberTimingCorrections.isValid());
    CSCDetId idraw = mapper_->rawCSCDetId(id);
    int index = indexer_->chamberIndex(idraw) - 1;  // NOTE THE MINUS ONE!
    return float(theChamberTimingCorrections->item(index).anode_bx_offset * 1. /
                 theChamberTimingCorrections->precision());
  } else
    return 0;
}

float CSCConditions::averageGain() const {
  const float loEdge = 5.0;           // consider gains above this
  const float hiEdge = 10.0;          // consider gains below this
  const float loLimit = 6.0;          // lowest acceptable average gain
  const float hiLimit = 9.0;          // highest acceptable average gain
  const float expectedAverage = 7.5;  // default average gain

  if (theAverageGain > 0.)
    return theAverageGain;  // only recalculate if necessary

  int n_strip = 0;
  float gain_tot = 0.;

  CSCDBGains::GainContainer::const_iterator it;
  for (it = theGains->gains.begin(); it != theGains->gains.end(); ++it) {
    float the_gain = float(it->gain_slope) / theGains->scale();
    if (the_gain > loEdge && the_gain < hiEdge) {
      gain_tot += the_gain;
      ++n_strip;
    }
  }

  // Average gain
  if (n_strip > 0) {
    theAverageGain = gain_tot / n_strip;
  }

  // Average gain has been around 7.5 in real data
  if (theAverageGain < loLimit || theAverageGain > hiLimit) {
    //    LogTrace("CSC") << "Average CSC strip gain = "
    //                    << theAverageGain << "  is reset to expected value "
    //                    << expectedAverage;
    theAverageGain = expectedAverage;
  }

  return theAverageGain;
}
//
float CSCConditions::gasGainCorrection(const CSCDetId &id, int geomChannel, int iwiregroup) const {
  if (useGasGainCorrections()) {
    assert(theGasGainCorrections.isValid());
    CSCDetId idraw = mapper_->rawCSCDetId(id);
    int iraw = mapper_->rawStripChannel(id, geomChannel);
    int index = indexer_->gasGainIndex(idraw, iraw, iwiregroup) - 1;  // NOTE THE MINUS ONE!
    return float(theGasGainCorrections->value(index));
  } else {
    return 1.;
  }
}

int CSCConditions::channelFromStrip(const CSCDetId &id, int geomStrip) const {
  return mapper_->channelFromStrip(id, geomStrip);
}

int CSCConditions::rawStripChannel(const CSCDetId &id, int geomChannel) const {
  return mapper_->rawStripChannel(id, geomChannel);
}

void CSCConditions::print() const
//@@ HAS NOT BEEN UPDATED THROUGH SEVERAL VERSIONS OF THE CONDITIONS DATA
{
  /*
    std::cout << "SIZES: GAINS: " << theGains->gains.size()
              << "   PEDESTALS: " << thePedestals->pedestals.size()
              << "   NOISES "  << theNoiseMatrix->matrix.size() << std::endl;;

    std::map< int,std::vector<CSCDBGains::Item> >::const_iterator layerGainsItr
    = theGains->gains.begin(), lastGain = theGains->gains.end(); for( ;
    layerGainsItr != lastGain; ++layerGainsItr)
    {
      std::cout << "GAIN " << layerGainsItr->first
                << " STRIPS " << layerGainsItr->second.size() << " "
                << layerGainsItr->second[0].gain_slope
                << " " << layerGainsItr->second[0].gain_intercept << std::endl;
    }

    std::map< int,std::vector<CSCDBPedestals::Item> >::const_iterator
    pedestalItr = thePedestals->pedestals.begin(), lastPedestal =
    thePedestals->pedestals.end(); for( ; pedestalItr != lastPedestal;
    ++pedestalItr)
    {
      std::cout << "PEDS " << pedestalItr->first << " "
                << " STRIPS " << pedestalItr->second.size() << " ";
      for(int i = 1; i < 80; ++i)
      {
         std::cout << pedestalItr->second[i-1].rms << " " ;
       }
       std::cout << std::endl;
    }

    std::map< int,std::vector<CSCDBCrosstalk::Item> >::const_iterator
    crosstalkItr = theCrosstalk->crosstalk.begin(), lastCrosstalk =
    theCrosstalk->crosstalk.end(); for( ; crosstalkItr != lastCrosstalk;
    ++crosstalkItr)
    {
      std::cout << "XTALKS " << crosstalkItr->first
        << " STRIPS " << crosstalkItr->second.size() << " "
       << crosstalkItr->second[5].xtalk_slope_left << " "
       << crosstalkItr->second[5].xtalk_slope_right << " "
       << crosstalkItr->second[5].xtalk_intercept_left << " "
       << crosstalkItr->second[5].xtalk_intercept_right << std::endl;
    }
  */
}