CSCDigitizer/src/CSCBaseElectronicsSim.cc

0001 // This is CSCBaseElectronicsSim.cc
0002
0003 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0006 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0007 #include "SimMuon/CSCDigitizer/src/CSCBaseElectronicsSim.h"
0008 #include "SimMuon/CSCDigitizer/src/CSCDetectorHit.h"
0009
0010 #include "CLHEP/Random/RandGaussQ.h"
0011
0012 #include <algorithm>
0013 #include <list>
0014
0015 CSCBaseElectronicsSim::CSCBaseElectronicsSim(const edm::ParameterSet &p)
0016     : theSpecs(nullptr),
0017       theLayerGeometry(nullptr),
0018       theLayer(nullptr),
0019       theSignalMap(),
0020       theAmpResponse(),
0021       theBunchSpacing(25.),
0022       theNoiseWasAdded(false),
0023       nElements(0),
0024       theShapingTime(p.getParameter<int>("shapingTime")),
0025       thePeakTimeSigma(p.getParameter<double>("peakTimeSigma")),
0026       theBunchTimingOffsets(p.getParameter<std::vector<double>>("bunchTimingOffsets")),
0027       theSignalStartTime(p.getParameter<double>("signalStartTime")),
0028       theSignalStopTime(p.getParameter<double>("signalStopTime")),
0029       theSamplingTime(p.getParameter<double>("samplingTime")),
0030       theNumberOfSamples(static_cast<int>((theSignalStopTime - theSignalStartTime) / theSamplingTime)),
0031       theOffsetOfBxZero(p.getParameter<int>("timeBitForBxZero")),
0032       theSignalPropagationSpeed(p.getParameter<std::vector<double>>("signalSpeed")),
0033       theTimingCalibrationError(p.getParameter<std::vector<double>>("timingCalibrationError")),
0034       doNoise_(p.getParameter<bool>("doNoise")) {
0035   assert(theBunchTimingOffsets.size() == 11);
0036 }
0037
0038 CSCBaseElectronicsSim::~CSCBaseElectronicsSim() {}
0039
0040 void CSCBaseElectronicsSim::simulate(const CSCLayer *layer,
0041                                      const std::vector<CSCDetectorHit> &detectorHits,
0042                                      CLHEP::HepRandomEngine *engine) {
0043   theNoiseWasAdded = false;
0044
0045   {
0046     theSignalMap.clear();
0047     theDetectorHitMap.clear();
0048     setLayer(layer);
0049     // can we swap for efficiency?
0050     theDigiSimLinks = DigiSimLinks(layerId().rawId());
0051   }
0052
0053   {
0054     size_t nHits = detectorHits.size();
0055     // turn each detector hit into an analog signal
0056     for (size_t i = 0; i < nHits; ++i) {
0057       int element = readoutElement(detectorHits[i].getElement());
0058
0059       // skip if  hit element is not part of a readout element
0060       // e.g. wire in non-readout group
0061       if (element != 0)
0062         add(amplifySignal(detectorHits[i]), engine);
0063     }
0064   }
0065
0066   {
0067     if (doNoise_) {
0068       addNoise(engine);
0069     }
0070   }
0071 }
0072
0073 void CSCBaseElectronicsSim::setLayer(const CSCLayer *layer) {
0074   // fill the specs member data
0075   theSpecs = layer->chamber()->specs();
0076   theLayerGeometry = layer->geometry();
0077
0078   theLayer = layer;
0079   theLayerId = CSCDetId(theLayer->geographicalId().rawId());
0080   initParameters();
0081 }
0082
0083 void CSCBaseElectronicsSim::fillAmpResponse() {
0084   std::vector<float> ampBinValues(theNumberOfSamples);
0085   int i = 0;
0086   for (; i < theNumberOfSamples; ++i) {
0087     ampBinValues[i] = calculateAmpResponse(i * theSamplingTime);
0088     // truncate any entries that are trivially small
0089     if (i > 5 && ampBinValues[i] < 0.000001)
0090       break;
0091   }
0092   ampBinValues.resize(i);
0093   theAmpResponse = CSCAnalogSignal(0, theSamplingTime, ampBinValues, 1., 0.);
0094
0095   LogTrace("CSCBaseElectronicsSim") << "CSCBaseElectronicsSim: dump of theAmpResponse follows...\n" << theAmpResponse;
0096 }
0097
0098 CSCAnalogSignal CSCBaseElectronicsSim::amplifySignal(const CSCDetectorHit &detectorHit) {
0099   int element = readoutElement(detectorHit.getElement());
0100
0101   float readoutTime = detectorHit.getTime() + signalDelay(element, detectorHit.getPosition());
0102
0103   // start from the amp response, and modify it.
0104   CSCAnalogSignal thisSignal(theAmpResponse);
0105   thisSignal *= detectorHit.getCharge();
0106   thisSignal.setTimeOffset(readoutTime);
0107   thisSignal.setElement(element);
0108   // keep track of links between digis and hits
0109   theDetectorHitMap.insert(DetectorHitMap::value_type(channelIndex(element), detectorHit));
0110   return thisSignal;
0111 }
0112
0113 CSCAnalogSignal CSCBaseElectronicsSim::makeNoiseSignal(int element, CLHEP::HepRandomEngine *) {
0114   std::vector<float> binValues(theNumberOfSamples);
0115   // default is empty
0116   return CSCAnalogSignal(element, theSamplingTime, binValues, 0., theSignalStartTime);
0117 }
0118
0119 void CSCBaseElectronicsSim::addNoise(CLHEP::HepRandomEngine *engine) {
0120   for (CSCSignalMap::iterator mapI = theSignalMap.begin(); mapI != theSignalMap.end(); ++mapI) {
0121     // superimpose electronics noise
0122     (*mapI).second.superimpose(makeNoiseSignal((*mapI).first, engine));
0123     // DON'T do amp gain variations.  Handled in strips by calibration code
0124     // and variations in the shaper peaking time.
0125     double timeOffset = CLHEP::RandGaussQ::shoot(engine, (*mapI).second.getTimeOffset(), thePeakTimeSigma);
0126     (*mapI).second.setTimeOffset(timeOffset);
0127   }
0128   theNoiseWasAdded = true;
0129 }
0130
0131 CSCAnalogSignal &CSCBaseElectronicsSim::find(int element, CLHEP::HepRandomEngine *engine) {
0132   if (element <= 0 || element > nElements) {
0133     LogTrace("CSCBaseElectronicsSim") << "CSCBaseElectronicsSim: bad element = " << element << ". There are "
0134                                       << nElements << " elements.";
0135     edm::LogError("Error in CSCBaseElectronicsSim:  element out of bounds");
0136   }
0137   CSCSignalMap::iterator signalMapItr = theSignalMap.find(element);
0138   if (signalMapItr == theSignalMap.end()) {
0139     CSCAnalogSignal newSignal;
0140     if (theNoiseWasAdded) {
0141       newSignal = makeNoiseSignal(element, engine);
0142     } else {
0143       std::vector<float> emptyV(theNumberOfSamples);
0144       newSignal = CSCAnalogSignal(element, theSamplingTime, emptyV, 0., theSignalStartTime);
0145     }
0146     signalMapItr = theSignalMap.insert(std::pair<int, CSCAnalogSignal>(element, newSignal)).first;
0147   }
0148   return (*signalMapItr).second;
0149 }
0150
0151 CSCAnalogSignal &CSCBaseElectronicsSim::add(const CSCAnalogSignal &signal, CLHEP::HepRandomEngine *engine) {
0152   int element = signal.getElement();
0153   CSCAnalogSignal &newSignal = find(element, engine);
0154   newSignal.superimpose(signal);
0155   return newSignal;
0156 }
0157
0158 float CSCBaseElectronicsSim::signalDelay(int element, float pos) const {
0159   // readout is on top edge of chamber for strips, right edge
0160   // for wires.
0161   // zero calibrated to chamber center
0162   float distance = -1. * pos;
0163   float speed = theSignalPropagationSpeed[theSpecs->chamberType()];
0164   return distance / speed;
0165 }
0166
0167 void CSCBaseElectronicsSim::addLinks(int channelIndex) {
0168   std::pair<DetectorHitMap::iterator, DetectorHitMap::iterator> channelHitItr =
0169       theDetectorHitMap.equal_range(channelIndex);
0170
0171   // find the fraction contribution for each SimTrack
0172   std::map<int, float> simTrackChargeMap;
0173   std::map<int, EncodedEventId> eventIdMap;
0174   float totalCharge = 0;
0175   for (DetectorHitMap::iterator hitItr = channelHitItr.first; hitItr != channelHitItr.second; ++hitItr) {
0176     const PSimHit *hit = hitItr->second.getSimHit();
0177     // might be zero for unit tests and such
0178     if (hit != nullptr) {
0179       int simTrackId = hitItr->second.getSimHit()->trackId();
0180       float charge = hitItr->second.getCharge();
0181       std::map<int, float>::iterator chargeItr = simTrackChargeMap.find(simTrackId);
0182       if (chargeItr == simTrackChargeMap.end()) {
0183         simTrackChargeMap[simTrackId] = charge;
0184         eventIdMap[simTrackId] = hit->eventId();
0185       } else {
0186         chargeItr->second += charge;
0187       }
0188       totalCharge += charge;
0189     }
0190   }
0191
0192   for (std::map<int, float>::iterator chargeItr = simTrackChargeMap.begin(); chargeItr != simTrackChargeMap.end();
0193        ++chargeItr) {
0194     int simTrackId = chargeItr->first;
0195     theDigiSimLinks.push_back(
0196         StripDigiSimLink(channelIndex, simTrackId, eventIdMap[simTrackId], chargeItr->second / totalCharge));
0197   }
0198 }