CSCBaseElectronicsSim.cc

Go to the documentation of this file.

// This is CSCBaseElectronicsSim.cc

#include "CLHEP/Units/GlobalPhysicalConstants.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/CSCGeometry/interface/CSCLayer.h"
#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
#include "SimMuon/CSCDigitizer/src/CSCBaseElectronicsSim.h"
#include "SimMuon/CSCDigitizer/src/CSCDetectorHit.h"

#include "CLHEP/Random/RandGaussQ.h"

#include <algorithm>
#include <list>

CSCBaseElectronicsSim::CSCBaseElectronicsSim(const edm::ParameterSet &p)
    : theSpecs(nullptr),
      theLayerGeometry(nullptr),
      theLayer(nullptr),
      theSignalMap(),
      theAmpResponse(),
      theBunchSpacing(25.),
      theNoiseWasAdded(false),
      nElements(0),
      theShapingTime(p.getParameter<int>("shapingTime")),
      thePeakTimeSigma(p.getParameter<double>("peakTimeSigma")),
      theBunchTimingOffsets(p.getParameter<std::vector<double>>("bunchTimingOffsets")),
      theSignalStartTime(p.getParameter<double>("signalStartTime")),
      theSignalStopTime(p.getParameter<double>("signalStopTime")),
      theSamplingTime(p.getParameter<double>("samplingTime")),
      theNumberOfSamples(static_cast<int>((theSignalStopTime - theSignalStartTime) / theSamplingTime)),
      theOffsetOfBxZero(p.getParameter<int>("timeBitForBxZero")),
      theSignalPropagationSpeed(p.getParameter<std::vector<double>>("signalSpeed")),
      theTimingCalibrationError(p.getParameter<std::vector<double>>("timingCalibrationError")),
      doNoise_(p.getParameter<bool>("doNoise")) {
  assert(theBunchTimingOffsets.size() == 11);
}

CSCBaseElectronicsSim::~CSCBaseElectronicsSim() {}

void CSCBaseElectronicsSim::simulate(const CSCLayer *layer,
                                     const std::vector<CSCDetectorHit> &detectorHits,
                                     CLHEP::HepRandomEngine *engine) {
  theNoiseWasAdded = false;

  {
    theSignalMap.clear();
    theDetectorHitMap.clear();
    setLayer(layer);
    // can we swap for efficiency?
    theDigiSimLinks = DigiSimLinks(layerId().rawId());
  }

  {
    size_t nHits = detectorHits.size();
    // turn each detector hit into an analog signal
    for (size_t i = 0; i < nHits; ++i) {
      int element = readoutElement(detectorHits[i].getElement());

      // skip if  hit element is not part of a readout element
      // e.g. wire in non-readout group
      if (element != 0)
        add(amplifySignal(detectorHits[i]), engine);
    }
  }

  {
    if (doNoise_) {
      addNoise(engine);
    }
  }
}

void CSCBaseElectronicsSim::setLayer(const CSCLayer *layer) {
  // fill the specs member data
  theSpecs = layer->chamber()->specs();
  theLayerGeometry = layer->geometry();

  theLayer = layer;
  theLayerId = CSCDetId(theLayer->geographicalId().rawId());
  initParameters();
}

void CSCBaseElectronicsSim::fillAmpResponse() {
  std::vector<float> ampBinValues(theNumberOfSamples);
  int i = 0;
  for (; i < theNumberOfSamples; ++i) {
    ampBinValues[i] = calculateAmpResponse(i * theSamplingTime);
    // truncate any entries that are trivially small
    if (i > 5 && ampBinValues[i] < 0.000001)
      break;
  }
  ampBinValues.resize(i);
  theAmpResponse = CSCAnalogSignal(0, theSamplingTime, ampBinValues, 1., 0.);

  LogTrace("CSCBaseElectronicsSim") << "CSCBaseElectronicsSim: dump of theAmpResponse follows...\n" << theAmpResponse;
}

CSCAnalogSignal CSCBaseElectronicsSim::amplifySignal(const CSCDetectorHit &detectorHit) {
  int element = readoutElement(detectorHit.getElement());

  float readoutTime = detectorHit.getTime() + signalDelay(element, detectorHit.getPosition());

  // start from the amp response, and modify it.
  CSCAnalogSignal thisSignal(theAmpResponse);
  thisSignal *= detectorHit.getCharge();
  thisSignal.setTimeOffset(readoutTime);
  thisSignal.setElement(element);
  // keep track of links between digis and hits
  theDetectorHitMap.insert(DetectorHitMap::value_type(channelIndex(element), detectorHit));
  return thisSignal;
}

CSCAnalogSignal CSCBaseElectronicsSim::makeNoiseSignal(int element, CLHEP::HepRandomEngine *) {
  std::vector<float> binValues(theNumberOfSamples);
  // default is empty
  return CSCAnalogSignal(element, theSamplingTime, binValues, 0., theSignalStartTime);
}

void CSCBaseElectronicsSim::addNoise(CLHEP::HepRandomEngine *engine) {
  for (CSCSignalMap::iterator mapI = theSignalMap.begin(); mapI != theSignalMap.end(); ++mapI) {
    // superimpose electronics noise
    (*mapI).second.superimpose(makeNoiseSignal((*mapI).first, engine));
    // DON'T do amp gain variations.  Handled in strips by calibration code
    // and variations in the shaper peaking time.
    double timeOffset = CLHEP::RandGaussQ::shoot(engine, (*mapI).second.getTimeOffset(), thePeakTimeSigma);
    (*mapI).second.setTimeOffset(timeOffset);
  }
  theNoiseWasAdded = true;
}

CSCAnalogSignal &CSCBaseElectronicsSim::find(int element, CLHEP::HepRandomEngine *engine) {
  if (element <= 0 || element > nElements) {
    LogTrace("CSCBaseElectronicsSim") << "CSCBaseElectronicsSim: bad element = " << element << ". There are "
                                      << nElements << " elements.";
    edm::LogError("Error in CSCBaseElectronicsSim:  element out of bounds");
  }
  CSCSignalMap::iterator signalMapItr = theSignalMap.find(element);
  if (signalMapItr == theSignalMap.end()) {
    CSCAnalogSignal newSignal;
    if (theNoiseWasAdded) {
      newSignal = makeNoiseSignal(element, engine);
    } else {
      std::vector<float> emptyV(theNumberOfSamples);
      newSignal = CSCAnalogSignal(element, theSamplingTime, emptyV, 0., theSignalStartTime);
    }
    signalMapItr = theSignalMap.insert(std::pair<int, CSCAnalogSignal>(element, newSignal)).first;
  }
  return (*signalMapItr).second;
}

CSCAnalogSignal &CSCBaseElectronicsSim::add(const CSCAnalogSignal &signal, CLHEP::HepRandomEngine *engine) {
  int element = signal.getElement();
  CSCAnalogSignal &newSignal = find(element, engine);
  newSignal.superimpose(signal);
  return newSignal;
}

float CSCBaseElectronicsSim::signalDelay(int element, float pos) const {
  // readout is on top edge of chamber for strips, right edge
  // for wires.
  // zero calibrated to chamber center
  float distance = -1. * pos;
  float speed = theSignalPropagationSpeed[theSpecs->chamberType()];
  return distance / speed;
}

void CSCBaseElectronicsSim::addLinks(int channelIndex) {
  std::pair<DetectorHitMap::iterator, DetectorHitMap::iterator> channelHitItr =
      theDetectorHitMap.equal_range(channelIndex);

  // find the fraction contribution for each SimTrack
  std::map<int, float> simTrackChargeMap;
  std::map<int, EncodedEventId> eventIdMap;
  float totalCharge = 0;
  for (DetectorHitMap::iterator hitItr = channelHitItr.first; hitItr != channelHitItr.second; ++hitItr) {
    const PSimHit *hit = hitItr->second.getSimHit();
    // might be zero for unit tests and such
    if (hit != nullptr) {
      int simTrackId = hitItr->second.getSimHit()->trackId();
      float charge = hitItr->second.getCharge();
      std::map<int, float>::iterator chargeItr = simTrackChargeMap.find(simTrackId);
      if (chargeItr == simTrackChargeMap.end()) {
        simTrackChargeMap[simTrackId] = charge;
        eventIdMap[simTrackId] = hit->eventId();
      } else {
        chargeItr->second += charge;
      }
      totalCharge += charge;
    }
  }

  for (std::map<int, float>::iterator chargeItr = simTrackChargeMap.begin(); chargeItr != simTrackChargeMap.end();
       ++chargeItr) {
    int simTrackId = chargeItr->first;
    theDigiSimLinks.push_back(
        StripDigiSimLink(channelIndex, simTrackId, eventIdMap[simTrackId], chargeItr->second / totalCharge));
  }
}