Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:34:55

0001 
0002 #include "CalibMuon/CSCCalibration/interface/CSCConditions.h"
0003 
0004 #include "CalibMuon/CSCCalibration/interface/CSCChannelMapperBase.h"
0005 #include "CalibMuon/CSCCalibration/interface/CSCIndexerBase.h"
0006 
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/ESWatcher.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 
0014 #include "CondFormats/CSCObjects/interface/CSCChamberTimeCorrections.h"
0015 #include "CondFormats/CSCObjects/interface/CSCDBChipSpeedCorrection.h"
0016 #include "CondFormats/CSCObjects/interface/CSCDBCrosstalk.h"
0017 #include "CondFormats/CSCObjects/interface/CSCDBGains.h"
0018 #include "CondFormats/CSCObjects/interface/CSCDBGasGainCorrection.h"
0019 #include "CondFormats/CSCObjects/interface/CSCDBPedestals.h"
0020 
0021 #include "CondFormats/CSCObjects/interface/CSCBadChambers.h"
0022 #include "CondFormats/CSCObjects/interface/CSCBadStrips.h"
0023 #include "CondFormats/CSCObjects/interface/CSCBadWires.h"
0024 
0025 CSCConditions::CSCConditions(const edm::ParameterSet &ps, edm::ConsumesCollector cc)
0026     : theGains(),
0027       theCrosstalk(),
0028       thePedestals(),
0029       theNoiseMatrix(),
0030       theBadStrips(),
0031       theBadWires(),
0032       theBadChambers(),
0033       theChipCorrections(),
0034       theChamberTimingCorrections(),
0035       theGasGainCorrections(),
0036       indexer_(nullptr),
0037       mapper_(nullptr),
0038       readBadChannels_(false),
0039       readBadChambers_(false),
0040       useTimingCorrections_(false),
0041       useGasGainCorrections_(false),
0042       idOfBadChannelWords_(CSCDetId()),
0043       badStripWord_(0),
0044       badWireWord_(0),
0045       theAverageGain(-1.0) {
0046   readBadChannels_ = ps.getParameter<bool>("readBadChannels");
0047   readBadChambers_ = ps.getParameter<bool>("readBadChambers");
0048   useTimingCorrections_ = ps.getParameter<bool>("CSCUseTimingCorrections");
0049   useGasGainCorrections_ = ps.getParameter<bool>("CSCUseGasGainCorrections");
0050   indexerToken_ = cc.esConsumes<CSCIndexerBase, CSCIndexerRecord>();
0051   mapperToken_ = cc.esConsumes<CSCChannelMapperBase, CSCChannelMapperRecord>();
0052   gainsToken_ = cc.esConsumes<CSCDBGains, CSCDBGainsRcd>();
0053   crosstalkToken_ = cc.esConsumes<CSCDBCrosstalk, CSCDBCrosstalkRcd>();
0054   pedestalsToken_ = cc.esConsumes<CSCDBPedestals, CSCDBPedestalsRcd>();
0055   noiseMatrixToken_ = cc.esConsumes<CSCDBNoiseMatrix, CSCDBNoiseMatrixRcd>();
0056   if (useTimingCorrections()) {
0057     chipCorrectionsToken_ = cc.esConsumes<CSCDBChipSpeedCorrection, CSCDBChipSpeedCorrectionRcd>();
0058     chamberTimingCorrectionsToken_ = cc.esConsumes<CSCChamberTimeCorrections, CSCChamberTimeCorrectionsRcd>();
0059   }
0060   if (readBadChannels()) {
0061     badStripsToken_ = cc.esConsumes<CSCBadStrips, CSCBadStripsRcd>();
0062     badWiresToken_ = cc.esConsumes<CSCBadWires, CSCBadWiresRcd>();
0063   }
0064   if (readBadChambers()) {
0065     badChambersToken_ = cc.esConsumes<CSCBadChambers, CSCBadChambersRcd>();
0066   }
0067   if (useGasGainCorrections()) {
0068     gasGainCorrectionsToken_ = cc.esConsumes<CSCDBGasGainCorrection, CSCDBGasGainCorrectionRcd>();
0069   }
0070 }
0071 
0072 CSCConditions::~CSCConditions() {}
0073 
0074 void CSCConditions::initializeEvent(const edm::EventSetup &es) {
0075   // Algorithms
0076   indexer_ = es.getHandle(indexerToken_);
0077   mapper_ = es.getHandle(mapperToken_);
0078 
0079   // Strip gains
0080   theGains = es.getHandle(gainsToken_);
0081   // Strip X-talk
0082   theCrosstalk = es.getHandle(crosstalkToken_);
0083   // Strip pedestals
0084   thePedestals = es.getHandle(pedestalsToken_);
0085   // Strip autocorrelation noise matrix
0086   theNoiseMatrix = es.getHandle(noiseMatrixToken_);
0087 
0088   if (useTimingCorrections()) {
0089     // Buckeye chip speeds
0090     theChipCorrections = es.getHandle(chipCorrectionsToken_);
0091     // Cable lengths from chambers to peripheral crate and additional chamber
0092     // level timing correction
0093     theChamberTimingCorrections = es.getHandle(chamberTimingCorrectionsToken_);
0094   }
0095 
0096   if (readBadChannels()) {
0097     // Bad strip channels
0098     theBadStrips = es.getHandle(badStripsToken_);
0099     // Bad wiregroup channels
0100     theBadWires = es.getHandle(badWiresToken_);
0101 
0102     //@@    if( badStripsWatcher_.check( es ) ) {
0103     //      fillBadStripWords();
0104     //@@    }
0105     //@@    if( badWiresWatcher_.check( es ) ) {
0106     //      fillBadWireWords();
0107     //@    }
0108   }
0109 
0110   // Has GainsRcd changed?
0111   if (gainsWatcher_.check(es)) {  // Yes...
0112     theAverageGain = -1.0;        // ...reset, so next access will recalculate it
0113   }
0114 
0115   if (readBadChambers()) {
0116     // Entire bad chambers
0117     theBadChambers = es.getHandle(badChambersToken_);
0118   }
0119 
0120   if (useGasGainCorrections()) {
0121     theGasGainCorrections = es.getHandle(gasGainCorrectionsToken_);
0122   }
0123 
0124   //  print();
0125 }
0126 
0127 void CSCConditions::fillBadChannelWords(const CSCDetId &id) {
0128   // input CSCDetId is expected to be an offline value i.e. different for ME1/1A
0129   // and ME1/1B
0130 
0131   // Only update content if necessary
0132   if (id != idOfBadChannelWords()) {
0133     // store offline CSCDetId for the two bad channel words
0134     setIdOfBadChannelWords(id);
0135 
0136     // reset to all zeroes
0137     badStripWord_.reset();
0138     badWireWord_.reset();
0139 
0140     if (readBadChannels()) {
0141       // convert to online CSCDetId since that is how conditions data are stored
0142       CSCDetId idraw = mapper_->rawCSCDetId(id);
0143       fillBadStripWord(idraw);
0144       fillBadWireWord(idraw);
0145     }
0146   }
0147 }
0148 
0149 /// Next function private
0150 
0151 void CSCConditions::fillBadStripWord(const CSCDetId &id) {
0152   // Input CSCDetId is expected to be a 'raw' value
0153 
0154   // Find linear index of chamber for input CSCDetId
0155   int inputIndex = indexer_->chamberIndex(id);
0156   short inputLayer = id.layer();
0157 
0158   // Does this chamber occur in bad channel list? If so, unpack its bad channels
0159 
0160   // chambers is a vector<BadChamber>
0161   // channels is a vector<BadChannel>
0162   // Each BadChamber contains its index (1-468 or 540 w. ME42), the no. of bad
0163   // channels, and the index within vector<BadChannel> where this chamber's bad
0164   // channels start.
0165 
0166   for (size_t i = 0; i < theBadStrips->chambers.size(); ++i) {  // loop over bad chambers
0167     int indexc = theBadStrips->chambers[i].chamber_index;
0168     if (indexc != inputIndex)
0169       continue;  // next iteration if not a match
0170 
0171     int start = theBadStrips->chambers[i].pointer;
0172     int nbad = theBadStrips->chambers[i].bad_channels;
0173 
0174     for (int j = start - 1; j < start - 1 + nbad; ++j) {  // bad channels in this chamber
0175       short lay = theBadStrips->channels[j].layer;        // value 1-6
0176       if (lay != inputLayer)
0177         continue;
0178 
0179       short chan = theBadStrips->channels[j].channel;  // value 1-80 (->112 for unganged ME1/1A)
0180       // Flags so far unused (and unset in conditins data)
0181       //    short f1 = theBadStrips->channels[j].flag1;
0182       //    short f2 = theBadStrips->channels[j].flag2;
0183       //    short f3 = theBadStrips->channels[j].flag3;
0184       badStripWord_.set(chan - 1, true);  // set bit 0-79 (111) in 80 (112)-bit
0185                                           // bitset representing this layer
0186     }  // j
0187   }  // i
0188 }
0189 
0190 void CSCConditions::fillBadWireWord(const CSCDetId &id) {
0191   // Input CSCDetId is expected to be a 'raw' value
0192 
0193   // Find linear index of chamber for input CSCDetId
0194   int inputIndex = indexer_->chamberIndex(id);
0195   short inputLayer = id.layer();
0196 
0197   // unpack what we've read from theBadWires
0198 
0199   for (size_t i = 0; i < theBadWires->chambers.size(); ++i) {  // loop over bad chambers
0200     int indexc = theBadWires->chambers[i].chamber_index;
0201 
0202     if (indexc != inputIndex)
0203       continue;  // next iteration if not a match
0204 
0205     int start = theBadWires->chambers[i].pointer;
0206     int nbad = theBadWires->chambers[i].bad_channels;
0207 
0208     for (int j = start - 1; j < start - 1 + nbad; ++j) {  // bad channels in this chamber
0209       short lay = theBadWires->channels[j].layer;         // value 1-6
0210       if (lay != inputLayer)
0211         continue;
0212 
0213       short chan = theBadWires->channels[j].channel;  // value 1-112
0214       //    short f1 = theBadWires->channels[j].flag1;
0215       //    short f2 = theBadWires->channels[j].flag2;
0216       //    short f3 = theBadWires->channels[j].flag3;
0217       badWireWord_.set(chan - 1,
0218                        true);  // set bit 0-111 in 112-bit bitset representing this layer
0219     }  // j
0220   }  // i
0221 }
0222 
0223 bool CSCConditions::isInBadChamber(const CSCDetId &id) const {
0224   //@@ We do not consider the possibility of having ME1/1A & ME1/1B
0225   // independently 'bad'.
0226   //@@ To do that we would need to define separate chamber indexes for ME1/1A &
0227   // ME1/1B.
0228 
0229   if (readBadChambers()) {
0230     CSCDetId idraw = mapper_->rawCSCDetId(id);
0231     int index = indexer_->chamberIndex(idraw);
0232     return theBadChambers->isInBadChamber(index);
0233   } else
0234     return false;
0235 }
0236 
0237 float CSCConditions::gain(const CSCDetId &id, int geomChannel) const {
0238   assert(theGains.isValid());
0239   CSCDetId idraw = mapper_->rawCSCDetId(id);
0240   int iraw = mapper_->rawStripChannel(id, geomChannel);
0241   int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
0242   return float(theGains->gain(index)) / theGains->scale();
0243 }
0244 
0245 float CSCConditions::pedestal(const CSCDetId &id, int geomChannel) const {
0246   assert(thePedestals.isValid());
0247   CSCDetId idraw = mapper_->rawCSCDetId(id);
0248   int iraw = mapper_->rawStripChannel(id, geomChannel);
0249   int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
0250   return float(thePedestals->pedestal(index)) / thePedestals->scale_ped();
0251 }
0252 
0253 float CSCConditions::pedestalSigma(const CSCDetId &id, int geomChannel) const {
0254   assert(thePedestals.isValid());
0255   CSCDetId idraw = mapper_->rawCSCDetId(id);
0256   int iraw = mapper_->rawStripChannel(id, geomChannel);
0257   int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
0258   return float(thePedestals->pedestal_rms(index)) / thePedestals->scale_rms();
0259 }
0260 
0261 float CSCConditions::crosstalkIntercept(const CSCDetId &id, int geomChannel, bool leftRight) const {
0262   assert(theCrosstalk.isValid());
0263   CSCDetId idraw = mapper_->rawCSCDetId(id);
0264   int iraw = mapper_->rawStripChannel(id, geomChannel);
0265   int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
0266   // resistive fraction is at the peak, where t=0
0267   return leftRight ? float(theCrosstalk->rinter(index)) / theCrosstalk->iscale()
0268                    : float(theCrosstalk->linter(index)) / theCrosstalk->iscale();
0269 }
0270 
0271 float CSCConditions::crosstalkSlope(const CSCDetId &id, int geomChannel, bool leftRight) const {
0272   assert(theCrosstalk.isValid());
0273   CSCDetId idraw = mapper_->rawCSCDetId(id);
0274   int iraw = mapper_->rawStripChannel(id, geomChannel);
0275   int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
0276   // resistive fraction is at the peak, where t=0
0277   return leftRight ? float(theCrosstalk->rslope(index)) / theCrosstalk->sscale()
0278                    : float(theCrosstalk->lslope(index)) / theCrosstalk->sscale();
0279 }
0280 
0281 const CSCDBNoiseMatrix::Item &CSCConditions::noiseMatrix(const CSCDetId &id, int geomChannel) const {
0282   //@@ BEWARE - THIS FUNCTION DOES NOT APPLY SCALE FACTOR USED IN PACKING VALUES
0283   // IN CONDITIONS DATA
0284   //@@ MAY BE AN ERROR? WHO WOULD WANT ACCESS WITHOUT IT?
0285 
0286   assert(theNoiseMatrix.isValid());
0287   CSCDetId idraw = mapper_->rawCSCDetId(id);
0288   int iraw = mapper_->rawStripChannel(id, geomChannel);
0289   int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
0290   return theNoiseMatrix->item(index);
0291 }
0292 
0293 void CSCConditions::noiseMatrixElements(const CSCDetId &id, int geomChannel, std::vector<float> &me) const {
0294   assert(me.size() > 11);
0295   const CSCDBNoiseMatrix::Item &item = noiseMatrix(id, geomChannel);  // i.e. the function above
0296   me[0] = float(item.elem33) / theNoiseMatrix->scale();
0297   me[1] = float(item.elem34) / theNoiseMatrix->scale();
0298   me[2] = float(item.elem35) / theNoiseMatrix->scale();
0299   me[3] = float(item.elem44) / theNoiseMatrix->scale();
0300   me[4] = float(item.elem45) / theNoiseMatrix->scale();
0301   me[5] = float(item.elem46) / theNoiseMatrix->scale();
0302   me[6] = float(item.elem55) / theNoiseMatrix->scale();
0303   me[7] = float(item.elem56) / theNoiseMatrix->scale();
0304   me[8] = float(item.elem57) / theNoiseMatrix->scale();
0305   me[9] = float(item.elem66) / theNoiseMatrix->scale();
0306   me[10] = float(item.elem67) / theNoiseMatrix->scale();
0307   me[11] = float(item.elem77) / theNoiseMatrix->scale();
0308 }
0309 
0310 void CSCConditions::crossTalk(const CSCDetId &id, int geomChannel, std::vector<float> &ct) const {
0311   assert(theCrosstalk.isValid());
0312   CSCDetId idraw = mapper_->rawCSCDetId(id);
0313   int iraw = mapper_->rawStripChannel(id, geomChannel);
0314   int index = indexer_->stripChannelIndex(idraw, iraw) - 1;  // NOTE THE MINUS ONE!
0315 
0316   ct[0] = float(theCrosstalk->lslope(index)) / theCrosstalk->sscale();
0317   ct[1] = float(theCrosstalk->linter(index)) / theCrosstalk->iscale();
0318   ct[2] = float(theCrosstalk->rslope(index)) / theCrosstalk->sscale();
0319   ct[3] = float(theCrosstalk->rinter(index)) / theCrosstalk->iscale();
0320 }
0321 
0322 float CSCConditions::chipCorrection(const CSCDetId &id, int geomChannel) const {
0323   if (useTimingCorrections()) {
0324     assert(theChipCorrections.isValid());
0325     CSCDetId idraw = mapper_->rawCSCDetId(id);
0326     int iraw = mapper_->rawStripChannel(id, geomChannel);
0327     int ichip = indexer_->chipIndex(iraw);              // converts 1-80 to 1-5 (chip#, CFEB#)
0328     int index = indexer_->chipIndex(idraw, ichip) - 1;  // NOTE THE MINUS ONE!
0329     return float(theChipCorrections->value(index)) / theChipCorrections->scale();
0330   } else
0331     return 0;
0332 }
0333 float CSCConditions::chamberTimingCorrection(const CSCDetId &id) const {
0334   if (useTimingCorrections()) {
0335     assert(theChamberTimingCorrections.isValid());
0336     CSCDetId idraw = mapper_->rawCSCDetId(id);
0337     int index = indexer_->chamberIndex(idraw) - 1;  // NOTE THE MINUS ONE!
0338     return float(
0339         theChamberTimingCorrections->item(index).cfeb_tmb_skew_delay * 1. / theChamberTimingCorrections->precision() +
0340         theChamberTimingCorrections->item(index).cfeb_timing_corr * 1. / theChamberTimingCorrections->precision() +
0341         (theChamberTimingCorrections->item(index).cfeb_cable_delay * 25.));
0342   } else
0343     return 0;
0344 }
0345 float CSCConditions::anodeBXoffset(const CSCDetId &id) const {
0346   if (useTimingCorrections()) {
0347     assert(theChamberTimingCorrections.isValid());
0348     CSCDetId idraw = mapper_->rawCSCDetId(id);
0349     int index = indexer_->chamberIndex(idraw) - 1;  // NOTE THE MINUS ONE!
0350     return float(theChamberTimingCorrections->item(index).anode_bx_offset * 1. /
0351                  theChamberTimingCorrections->precision());
0352   } else
0353     return 0;
0354 }
0355 
0356 /// Return average strip gain for full CSC system. Lazy evaluation.
0357 /// Restrict averaging to gains between 5 and 10, and require average
0358 /// is between 6 or 9 otherwise fix it to 7.5.
0359 /// These values came from Dominique and Stan,
0360 float CSCConditions::averageGain() const {
0361   const float loEdge = 5.0;           // consider gains above this
0362   const float hiEdge = 10.0;          // consider gains below this
0363   const float loLimit = 6.0;          // lowest acceptable average gain
0364   const float hiLimit = 9.0;          // highest acceptable average gain
0365   const float expectedAverage = 7.5;  // default average gain
0366 
0367   if (theAverageGain > 0.)
0368     return theAverageGain;  // only recalculate if necessary
0369 
0370   int n_strip = 0;
0371   float gain_tot = 0.;
0372 
0373   CSCDBGains::GainContainer::const_iterator it;
0374   for (it = theGains->gains.begin(); it != theGains->gains.end(); ++it) {
0375     float the_gain = float(it->gain_slope) / theGains->scale();
0376     if (the_gain > loEdge && the_gain < hiEdge) {
0377       gain_tot += the_gain;
0378       ++n_strip;
0379     }
0380   }
0381 
0382   // Average gain
0383   if (n_strip > 0) {
0384     theAverageGain = gain_tot / n_strip;
0385   }
0386 
0387   // Average gain has been around 7.5 in real data
0388   if (theAverageGain < loLimit || theAverageGain > hiLimit) {
0389     //    LogTrace("CSC") << "Average CSC strip gain = "
0390     //                    << theAverageGain << "  is reset to expected value "
0391     //                    << expectedAverage;
0392     theAverageGain = expectedAverage;
0393   }
0394 
0395   return theAverageGain;
0396 }
0397 //
0398 float CSCConditions::gasGainCorrection(const CSCDetId &id, int geomChannel, int iwiregroup) const {
0399   if (useGasGainCorrections()) {
0400     assert(theGasGainCorrections.isValid());
0401     CSCDetId idraw = mapper_->rawCSCDetId(id);
0402     int iraw = mapper_->rawStripChannel(id, geomChannel);
0403     int index = indexer_->gasGainIndex(idraw, iraw, iwiregroup) - 1;  // NOTE THE MINUS ONE!
0404     return float(theGasGainCorrections->value(index));
0405   } else {
0406     return 1.;
0407   }
0408 }
0409 
0410 int CSCConditions::channelFromStrip(const CSCDetId &id, int geomStrip) const {
0411   return mapper_->channelFromStrip(id, geomStrip);
0412 }
0413 
0414 int CSCConditions::rawStripChannel(const CSCDetId &id, int geomChannel) const {
0415   return mapper_->rawStripChannel(id, geomChannel);
0416 }
0417 
0418 void CSCConditions::print() const
0419 //@@ HAS NOT BEEN UPDATED THROUGH SEVERAL VERSIONS OF THE CONDITIONS DATA
0420 {
0421   /*
0422     std::cout << "SIZES: GAINS: " << theGains->gains.size()
0423               << "   PEDESTALS: " << thePedestals->pedestals.size()
0424               << "   NOISES "  << theNoiseMatrix->matrix.size() << std::endl;;
0425 
0426     std::map< int,std::vector<CSCDBGains::Item> >::const_iterator layerGainsItr
0427     = theGains->gains.begin(), lastGain = theGains->gains.end(); for( ;
0428     layerGainsItr != lastGain; ++layerGainsItr)
0429     {
0430       std::cout << "GAIN " << layerGainsItr->first
0431                 << " STRIPS " << layerGainsItr->second.size() << " "
0432                 << layerGainsItr->second[0].gain_slope
0433                 << " " << layerGainsItr->second[0].gain_intercept << std::endl;
0434     }
0435 
0436     std::map< int,std::vector<CSCDBPedestals::Item> >::const_iterator
0437     pedestalItr = thePedestals->pedestals.begin(), lastPedestal =
0438     thePedestals->pedestals.end(); for( ; pedestalItr != lastPedestal;
0439     ++pedestalItr)
0440     {
0441       std::cout << "PEDS " << pedestalItr->first << " "
0442                 << " STRIPS " << pedestalItr->second.size() << " ";
0443       for(int i = 1; i < 80; ++i)
0444       {
0445          std::cout << pedestalItr->second[i-1].rms << " " ;
0446        }
0447        std::cout << std::endl;
0448     }
0449 
0450     std::map< int,std::vector<CSCDBCrosstalk::Item> >::const_iterator
0451     crosstalkItr = theCrosstalk->crosstalk.begin(), lastCrosstalk =
0452     theCrosstalk->crosstalk.end(); for( ; crosstalkItr != lastCrosstalk;
0453     ++crosstalkItr)
0454     {
0455       std::cout << "XTALKS " << crosstalkItr->first
0456         << " STRIPS " << crosstalkItr->second.size() << " "
0457        << crosstalkItr->second[5].xtalk_slope_left << " "
0458        << crosstalkItr->second[5].xtalk_slope_right << " "
0459        << crosstalkItr->second[5].xtalk_intercept_left << " "
0460        << crosstalkItr->second[5].xtalk_intercept_right << std::endl;
0461     }
0462   */
0463 }