CSCDigitizer/src/CSCStripElectronicsSim.cc

0001 #include "DataFormats/CSCDigi/interface/CSCComparatorDigi.h"
0002 #include "DataFormats/CSCDigi/interface/CSCStripDigi.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "Geometry/CSCGeometry/interface/CSCChamberSpecs.h"
0005 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0006 #include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
0007 #include "SimMuon/CSCDigitizer/src/CSCAnalogSignal.h"
0008 #include "SimMuon/CSCDigitizer/src/CSCCrosstalkGenerator.h"
0009 #include "SimMuon/CSCDigitizer/src/CSCDetectorHit.h"
0010 #include "SimMuon/CSCDigitizer/src/CSCStripConditions.h"
0011 #include "SimMuon/CSCDigitizer/src/CSCStripElectronicsSim.h"
0012
0013 #include <CLHEP/Random/RandGaussQ.h>
0014 #include <CLHEP/Units/GlobalPhysicalConstants.h>
0015 #include <CLHEP/Units/SystemOfUnits.h>
0016
0017 #include <boost/bind/bind.hpp>
0018 #include <cassert>
0019 #include <list>
0020
0021 using namespace boost::placeholders;
0022
0023 // This is CSCStripElectronicsSim.cc
0024
0025 CSCStripElectronicsSim::CSCStripElectronicsSim(const edm::ParameterSet &p)
0026     : CSCBaseElectronicsSim(p),
0027       theAmpResponse(theShapingTime, CSCStripAmpResponse::RADICAL),
0028       theComparatorThreshold(20.),
0029       theComparatorNoise(0.),
0030       theComparatorRMSOffset(2.),
0031       theComparatorSaturation(1057.),
0032       theComparatorWait(50.),
0033       theComparatorDeadTime(100.),
0034       theDaqDeadTime(200.),
0035       theTimingOffset(0.),
0036       nScaBins_(p.getParameter<int>("nScaBins")),
0037       doSuppression_(p.getParameter<bool>("doSuppression")),
0038       doCrosstalk_(p.getParameter<bool>("doCrosstalk")),
0039       theStripConditions(nullptr),
0040       theCrosstalkGenerator(nullptr),
0041       theComparatorClockJump(2),
0042       sca_time_bin_size(50.),
0043       sca_peak_bin(p.getParameter<int>("scaPeakBin")),
0044       theComparatorTimeBinOffset(p.getParameter<double>("comparatorTimeBinOffset")),
0045       theComparatorTimeOffset(p.getParameter<double>("comparatorTimeOffset")),
0046       theComparatorSamplingTime(p.getParameter<double>("comparatorSamplingTime")),
0047       theSCATimingOffsets(p.getParameter<std::vector<double>>("scaTimingOffsets")) {
0048   if (doCrosstalk_) {
0049     theCrosstalkGenerator = new CSCCrosstalkGenerator();
0050   }
0051
0052   fillAmpResponse();
0053 }
0054
0055 CSCStripElectronicsSim::~CSCStripElectronicsSim() {
0056   if (doCrosstalk_) {
0057     delete theCrosstalkGenerator;
0058   }
0059 }
0060
0061 void CSCStripElectronicsSim::initParameters() {
0062   nElements = theLayerGeometry->numberOfStrips();
0063   theComparatorThreshold = 20.;
0064   // selfTest();
0065
0066   // calculate the offset to the peak
0067   float averageDistance = theLayer->surface().position().mag();
0068   theAverageTimeOfFlight = averageDistance * CLHEP::cm / c_light;  // Units of c_light: mm/ns
0069   int chamberType = theSpecs->chamberType();
0070   theTimingOffset = theShapingTime + theAverageTimeOfFlight + theBunchTimingOffsets[chamberType];
0071   // TODO make sure config gets overridden
0072   theSignalStartTime = theTimingOffset - (sca_peak_bin - 1) * sca_time_bin_size;
0073   theSignalStopTime = theSignalStartTime + nScaBins_ * sca_time_bin_size;
0074   theNumberOfSamples = nScaBins_ * static_cast<int>(sca_time_bin_size / theSamplingTime);
0075 }
0076
0077 int CSCStripElectronicsSim::readoutElement(int strip) const { return theLayerGeometry->channel(strip); }
0078
0079 float CSCStripElectronicsSim::calculateAmpResponse(float t) const { return theAmpResponse.calculateAmpResponse(t); }
0080
0081 CSCAnalogSignal CSCStripElectronicsSim::makeNoiseSignal(int element, CLHEP::HepRandomEngine *engine) {
0082   std::vector<float> noiseBins(nScaBins_);
0083   CSCAnalogSignal tmpSignal(element, sca_time_bin_size, noiseBins);
0084   if (doNoise_) {
0085     theStripConditions->noisify(layerId(), tmpSignal, engine);
0086   }
0087   // now rebin it
0088   std::vector<float> binValues(theNumberOfSamples);
0089   for (int ibin = 0; ibin < theNumberOfSamples; ++ibin) {
0090     binValues[ibin] = tmpSignal.getValue(ibin * theSamplingTime);
0091   }
0092   CSCAnalogSignal finalSignal(element, theSamplingTime, binValues, 0., theSignalStartTime);
0093   return finalSignal;
0094 }
0095
0096 float CSCStripElectronicsSim::comparatorReading(const CSCAnalogSignal &signal,
0097                                                 float time,
0098                                                 CLHEP::HepRandomEngine *engine) const {
0099   return std::min(signal.getValue(time), theComparatorSaturation) +
0100          theComparatorRMSOffset * CLHEP::RandGaussQ::shoot(engine);
0101 }
0102
0103 void CSCStripElectronicsSim::runComparator(std::vector<CSCComparatorDigi> &result, CLHEP::HepRandomEngine *engine) {
0104   // first, make a list of all the comparators we actually
0105   // need to run
0106   std::list<int> comparatorsWithSignal;
0107   CSCSignalMap::iterator signalMapItr;
0108   for (signalMapItr = theSignalMap.begin(); signalMapItr != theSignalMap.end(); ++signalMapItr) {
0109     // Elements in signal map count from 1
0110     // 1,2->0,  3,4->1,  5,6->2, ...
0111     comparatorsWithSignal.push_back(((*signalMapItr).first - 1) / 2);
0112   }
0113   // no need to sort
0114   comparatorsWithSignal.unique();
0115   for (std::list<int>::iterator listItr = comparatorsWithSignal.begin(); listItr != comparatorsWithSignal.end();
0116        ++listItr) {
0117     int iComparator = *listItr;
0118     // find signal1 and signal2
0119     // iComparator counts from 0
0120     // icomp =0->1,2,  =1->3,4,  =2->5,6, ...
0121     const CSCAnalogSignal &signal1 = find(readoutElement(iComparator * 2 + 1), engine);
0122     const CSCAnalogSignal &signal2 = find(readoutElement(iComparator * 2 + 2), engine);
0123     for (float time = theSignalStartTime + theComparatorTimeOffset; time < theSignalStopTime - theComparatorWait;
0124          time += theComparatorSamplingTime) {
0125       if (comparatorReading(signal1, time, engine) > theComparatorThreshold ||
0126           comparatorReading(signal2, time, engine) > theComparatorThreshold) {
0127         // wait a bit, so we can run the comparator at the signal peak
0128         float comparatorTime = time;
0129         time += theComparatorWait;
0130
0131         float height1 = comparatorReading(signal1, time, engine);
0132         float height2 = comparatorReading(signal2, time, engine);
0133         int output = 0;
0134         int strip = 0;
0135         // distrip logic; comparator output is for pairs of strips:
0136         // hit  bin  dec
0137         // x--- 100   4
0138         // -x-- 101   5
0139         // --x- 110   6
0140         // ---x 111   7
0141         // just to prevent a copy
0142         const CSCAnalogSignal *mainSignal = nullptr;
0143         // pick the higher of the two strips in the pair
0144         if (height1 > height2) {
0145           mainSignal = &signal1;
0146           float leftStrip = 0.;
0147           if (iComparator > 0) {
0148             leftStrip = comparatorReading(find(readoutElement(iComparator * 2), engine), time, engine);
0149           }
0150           // if this strip is higher than either of its neighbors, make a
0151           // comparator digi
0152           if (leftStrip < height1 && height1 > theComparatorThreshold) {
0153             output = (leftStrip < height2);
0154             strip = iComparator * 2 + 1;
0155           }
0156         } else {
0157           mainSignal = &signal2;
0158           float rightStrip = 0.;
0159           if (iComparator * 2 + 3 <= nElements) {
0160             rightStrip = comparatorReading(find(readoutElement(iComparator * 2 + 3), engine), time, engine);
0161           }
0162           if (rightStrip < height2 && height2 > theComparatorThreshold) {
0163             output = (height1 < rightStrip);
0164             strip = iComparator * 2 + 2;
0165           }
0166         }
0167         if (strip != 0) {
0168           float bxFloat =
0169               (comparatorTime - theTimingOffset) / theBunchSpacing + theComparatorTimeBinOffset + theOffsetOfBxZero;
0170
0171           // Comparator digi as of Nov-2006 adapted to real data: time word has
0172           // 16 bits with set bit flagging appropriate bunch crossing, and bx 0
0173           // corresponding to 9th bit i.e.
0174
0175           //      1st bit set (bit 0) <-> bx -9
0176           //      2nd              1  <-> bx -8
0177           //      ...           ...       ....
0178           //      8th              9  <-> bx  0
0179           //      9th             10  <-> bx +1
0180           //      ...           ...       ....
0181           //     16th             15  <-> bx +6
0182
0183           // Parameter theOffsetOfBxZero = 9 @@WARNING! This offset may be
0184           // changed (hardware)!
0185
0186           int timeWord = 0;  // and this will remain if too early or late
0187           if ((bxFloat >= 0) && (bxFloat < 16))
0188             timeWord = (1 << static_cast<int>(bxFloat));  // set appropriate bit
0189
0190           CSCComparatorDigi newDigi(strip, output, timeWord);
0191           result.push_back(newDigi);
0192         }
0193
0194         // wait for the comparator to reset
0195         time += theComparatorDeadTime;
0196         // really should be zero, but strip signal doesn't go negative yet
0197         float resetThreshold = 1;
0198         while (time < theSignalStopTime && mainSignal->getValue(time) > resetThreshold) {
0199           time += theComparatorSamplingTime;
0200         }
0201
0202       }  // if over threshold
0203     }    // loop over time samples
0204   }      // loop over comparators
0205   // sort by time
0206   sort(result.begin(), result.end());
0207 }
0208
0209 std::list<int> CSCStripElectronicsSim::getKeyStrips(const std::vector<CSCComparatorDigi> &comparators) const {
0210   std::list<int> result;
0211   for (std::vector<CSCComparatorDigi>::const_iterator compItr = comparators.begin(); compItr != comparators.end();
0212        ++compItr) {
0213     if (std::abs(compItr->getTimeBin() - theOffsetOfBxZero) <= 2) {
0214       result.push_back(compItr->getStrip());
0215     }
0216   }
0217   // need sort for unique to work.
0218   result.sort();
0219   result.unique();
0220   return result;
0221 }
0222
0223 std::list<int> CSCStripElectronicsSim::getKeyStripsFromMC() const {
0224   // assumes the detector hit map is filled
0225   std::list<int> result;
0226   transform(theDetectorHitMap.begin(),
0227             theDetectorHitMap.end(),
0228             //   back_inserter(result),
0229             //   boost::bind(&DetectorHitMap::value_type::first,_1));
0230             // suggested code from Chris Jones
0231             back_inserter(result),
0232             std::bind(&DetectorHitMap::value_type::first, std::placeholders::_1));
0233   //  back_inserter(result), [](DetectorHitMap::value_type const& iValue) {
0234   //  return iValue.first; } );
0235   result.sort();
0236   result.unique();
0237   return result;
0238 }
0239
0240 std::list<int> CSCStripElectronicsSim::channelsToRead(const std::list<int> &keyStrips, int window) const {
0241   std::list<int> result;
0242   std::list<int>::const_iterator keyStripItr = keyStrips.begin();
0243   if (doSuppression_) {
0244     for (; keyStripItr != keyStrips.end(); ++keyStripItr) {
0245       // pick the five strips around the comparator
0246       for (int istrip = (*keyStripItr) - window; istrip <= (*keyStripItr) + window; ++istrip) {
0247         if (istrip > 0 && istrip <= nElements) {
0248           result.push_back(readoutElement(istrip));
0249         }
0250       }
0251     }
0252     result.sort();
0253     result.unique();
0254   } else {
0255     // read the whole CFEB, 16 strips
0256     std::list<int> cfebsToRead;
0257     for (; keyStripItr != keyStrips.end(); ++keyStripItr) {
0258       int cfeb = (readoutElement(*keyStripItr) - 1) / 16;
0259       cfebsToRead.push_back(cfeb);
0260       int remainder = (readoutElement(*keyStripItr) - 1) % 16;
0261       // if we're within 3 strips of an edge, take neighboring CFEB, too
0262       if (remainder < window && cfeb != 0) {
0263         cfebsToRead.push_back(cfeb - 1);
0264       }
0265       // the 'readouElement' makes it so that ME1/1 has just one CFEB
0266       int maxCFEBs = readoutElement(nElements) / 16 - 1;
0267       if (remainder >= 16 - window && cfeb != maxCFEBs) {
0268         cfebsToRead.push_back(cfeb + 1);
0269       }
0270     }
0271     cfebsToRead.sort();
0272     cfebsToRead.unique();
0273
0274     // now convert the CFEBS to strips
0275     for (std::list<int>::const_iterator cfebItr = cfebsToRead.begin(); cfebItr != cfebsToRead.end(); ++cfebItr) {
0276       for (int i = 1; i <= 16; ++i) {
0277         result.push_back((*cfebItr) * 16 + i);
0278       }
0279     }
0280   }
0281   return result;
0282 }
0283
0284 bool SortSignalsByTotal(const CSCAnalogSignal &s1, const CSCAnalogSignal &s2) {
0285   return (s1.getTotal() > s2.getTotal());
0286 }
0287
0288 void CSCStripElectronicsSim::fillDigis(CSCStripDigiCollection &digis,
0289                                        CSCComparatorDigiCollection &comparators,
0290                                        CLHEP::HepRandomEngine *engine) {
0291   if (doCrosstalk_) {
0292     addCrosstalk(engine);
0293   }
0294
0295   std::vector<CSCComparatorDigi> comparatorOutputs;
0296   runComparator(comparatorOutputs, engine);
0297   // copy these to the result
0298   if (!comparatorOutputs.empty()) {
0299     CSCComparatorDigiCollection::Range range(comparatorOutputs.begin(), comparatorOutputs.end());
0300     comparators.put(range, layerId());
0301   }
0302
0303   // std::list<int> keyStrips = getKeyStrips(comparatorOutputs);
0304   std::list<int> keyStrips = getKeyStripsFromMC();
0305   fillStripDigis(keyStrips, digis, engine);
0306 }
0307
0308 void CSCStripElectronicsSim::fillStripDigis(const std::list<int> &keyStrips,
0309                                             CSCStripDigiCollection &digis,
0310                                             CLHEP::HepRandomEngine *engine) {
0311   std::list<int> stripsToDo = channelsToRead(keyStrips, 3);
0312   std::vector<CSCStripDigi> stripDigis;
0313   stripDigis.reserve(stripsToDo.size());
0314   for (std::list<int>::const_iterator stripItr = stripsToDo.begin(); stripItr != stripsToDo.end(); ++stripItr) {
0315     createDigi(*stripItr, find(*stripItr, engine), stripDigis, engine);
0316   }
0317
0318   CSCStripDigiCollection::Range stripRange(stripDigis.begin(), stripDigis.end());
0319   digis.put(stripRange, layerId());
0320 }
0321
0322 void CSCStripElectronicsSim::addCrosstalk(CLHEP::HepRandomEngine *engine) {
0323   // this is needed so we can add a noise signal to the map
0324   // without messing up any iterators
0325   std::vector<CSCAnalogSignal> realSignals;
0326   realSignals.reserve(theSignalMap.size());
0327   CSCSignalMap::iterator mapI = theSignalMap.begin(), mapEnd = theSignalMap.end();
0328   for (; mapI != mapEnd; ++mapI) {
0329     realSignals.push_back((*mapI).second);
0330   }
0331   sort(realSignals.begin(), realSignals.end(), SortSignalsByTotal);
0332   std::vector<CSCAnalogSignal>::iterator realSignalItr = realSignals.begin(), realSignalsEnd = realSignals.end();
0333   for (; realSignalItr != realSignalsEnd; ++realSignalItr) {
0334     int thisStrip = (*realSignalItr).getElement();
0335     // add it to each neighbor
0336     if (thisStrip > 1) {
0337       int otherStrip = thisStrip - 1;
0338       addCrosstalk(*realSignalItr, thisStrip, otherStrip, engine);
0339     }
0340     if (thisStrip < nElements) {
0341       int otherStrip = thisStrip + 1;
0342       addCrosstalk(*realSignalItr, thisStrip, otherStrip, engine);
0343     }
0344   }
0345 }
0346
0347 void CSCStripElectronicsSim::addCrosstalk(const CSCAnalogSignal &signal,
0348                                           int thisStrip,
0349                                           int otherStrip,
0350                                           CLHEP::HepRandomEngine *engine) {
0351   float capacitiveCrosstalk, resistiveCrosstalk;
0352   bool leftRight = (otherStrip > thisStrip);
0353   theStripConditions->crosstalk(
0354       layerId(), thisStrip, theLayerGeometry->length(), leftRight, capacitiveCrosstalk, resistiveCrosstalk);
0355   theCrosstalkGenerator->setParameters(capacitiveCrosstalk, 0., resistiveCrosstalk);
0356   CSCAnalogSignal crosstalkSignal(theCrosstalkGenerator->getCrosstalk(signal));
0357   find(readoutElement(otherStrip), engine).superimpose(crosstalkSignal);
0358
0359   // Now subtract the crosstalk signal from the original signal
0360   crosstalkSignal *= -1.;
0361   find(thisStrip, engine).superimpose(crosstalkSignal);
0362 }
0363
0364 void CSCStripElectronicsSim::createDigi(int channel,
0365                                         const CSCAnalogSignal &signal,
0366                                         std::vector<CSCStripDigi> &result,
0367                                         CLHEP::HepRandomEngine *engine) {
0368   // fill in the sca information
0369   std::vector<int> scaCounts(nScaBins_);
0370
0371   float pedestal = theStripConditions->pedestal(layerId(), channel);
0372   float gain = theStripConditions->smearedGain(layerId(), channel, engine);
0373   int chamberType = theSpecs->chamberType();
0374   float timeSmearing = CLHEP::RandGaussQ::shoot(engine) * theTimingCalibrationError[chamberType];
0375   // undo the correction for TOF, instead, using some nominal
0376   // value from ME2/1
0377   float t0 = theSignalStartTime + theSCATimingOffsets[chamberType] + timeSmearing + 29. - theAverageTimeOfFlight;
0378   for (int scaBin = 0; scaBin < nScaBins_; ++scaBin) {
0379     float t = t0 + scaBin * sca_time_bin_size;
0380     scaCounts[scaBin] = static_cast<int>(pedestal + signal.getValue(t) * gain);
0381   }
0382   CSCStripDigi newDigi(channel, scaCounts);
0383
0384   // do saturation of 12-bit ADC
0385   doSaturation(newDigi);
0386
0387   result.push_back(newDigi);
0388   addLinks(channelIndex(channel));
0389   LogTrace("CSCStripElectronicsSim") << "CSCStripElectronicsSim: CSCStripDigi " << newDigi;
0390 }
0391
0392 void CSCStripElectronicsSim::doSaturation(CSCStripDigi &digi) {
0393   std::vector<int> scaCounts(digi.getADCCounts());
0394   for (unsigned scaBin = 0; scaBin < scaCounts.size(); ++scaBin) {
0395     scaCounts[scaBin] = std::min(scaCounts[scaBin], 4095);
0396   }
0397   digi.setADCCounts(scaCounts);
0398 }
0399
0400 void CSCStripElectronicsSim::fillMissingLayer(const CSCLayer *layer,
0401                                               const CSCComparatorDigiCollection &comparators,
0402                                               CSCStripDigiCollection &digis,
0403                                               CLHEP::HepRandomEngine *engine) {
0404   theSignalMap.clear();
0405   setLayer(layer);
0406   CSCDetId chamberId(theLayerId.chamberId());
0407   // find all comparator key strips in this chamber
0408   std::list<int> chamberKeyStrips;
0409   for (CSCComparatorDigiCollection::DigiRangeIterator comparatorItr = comparators.begin();
0410        comparatorItr != comparators.end();
0411        ++comparatorItr) {
0412     // could be more efficient
0413     if (CSCDetId((*comparatorItr).first).chamberId() == chamberId) {
0414       std::vector<CSCComparatorDigi> layerComparators((*comparatorItr).second.first, (*comparatorItr).second.second);
0415       std::list<int> layerKeyStrips = getKeyStrips(layerComparators);
0416       chamberKeyStrips.insert(chamberKeyStrips.end(), layerKeyStrips.begin(), layerKeyStrips.end());
0417     }
0418   }
0419   chamberKeyStrips.sort();
0420   chamberKeyStrips.unique();
0421   fillStripDigis(chamberKeyStrips, digis, engine);
0422 }
0423
0424 void CSCStripElectronicsSim::selfTest() const {
0425   // make sure the zero suppression algorithms work
0426   std::list<int> keyStrips, stripsRead;
0427   //
0428   bool isGanged = (readoutElement(nElements) == 16);
0429   keyStrips.push_back(readoutElement(19));
0430   keyStrips.push_back(readoutElement(30));
0431   keyStrips.push_back(readoutElement(32));
0432   stripsRead = channelsToRead(keyStrips, 3);
0433   if (doSuppression_) {
0434     unsigned int expectedSize = isGanged ? 10 : 12;
0435     assert(stripsRead.size() == expectedSize);
0436     assert(stripsRead.front() == readoutElement(17));
0437   } else {
0438     unsigned int expectedSize = isGanged ? 16 : 48;
0439     assert(stripsRead.size() == expectedSize);
0440     assert(stripsRead.front() == 1);
0441   }
0442 }