File indexing completed on 2023-03-17 11:25:18
0001
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
0050 theDigiSimLinks = DigiSimLinks(layerId().rawId());
0051 }
0052
0053 {
0054 size_t nHits = detectorHits.size();
0055
0056 for (size_t i = 0; i < nHits; ++i) {
0057 int element = readoutElement(detectorHits[i].getElement());
0058
0059
0060
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
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
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
0104 CSCAnalogSignal thisSignal(theAmpResponse);
0105 thisSignal *= detectorHit.getCharge();
0106 thisSignal.setTimeOffset(readoutTime);
0107 thisSignal.setElement(element);
0108
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
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
0122 (*mapI).second.superimpose(makeNoiseSignal((*mapI).first, engine));
0123
0124
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
0160
0161
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
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
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 }