File indexing completed on 2023-03-17 11:25:19
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0003 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0004 #include "SimMuon/CSCDigitizer/src/CSCDetectorHit.h"
0005 #include "SimMuon/CSCDigitizer/src/CSCDigitizer.h"
0006 #include "SimMuon/CSCDigitizer/src/CSCDriftSim.h"
0007 #include "SimMuon/CSCDigitizer/src/CSCNeutronReader.h"
0008 #include "SimMuon/CSCDigitizer/src/CSCStripConditions.h"
0009 #include "SimMuon/CSCDigitizer/src/CSCStripElectronicsSim.h"
0010 #include "SimMuon/CSCDigitizer/src/CSCStripHitSim.h"
0011 #include "SimMuon/CSCDigitizer/src/CSCWireElectronicsSim.h"
0012 #include "SimMuon/CSCDigitizer/src/CSCWireHitSim.h"
0013 #include <iostream>
0014
0015 CSCDigitizer::CSCDigitizer(const edm::ParameterSet &p)
0016 : theDriftSim(new CSCDriftSim()),
0017 theWireHitSim(new CSCWireHitSim(theDriftSim, p)),
0018 theStripHitSim(new CSCStripHitSim()),
0019 theWireElectronicsSim(new CSCWireElectronicsSim(p.getParameter<edm::ParameterSet>("wires"))),
0020 theStripElectronicsSim(new CSCStripElectronicsSim(p.getParameter<edm::ParameterSet>("strips"))),
0021 theNeutronReader(nullptr),
0022 theCSCGeometry(nullptr),
0023 theLayersNeeded(p.getParameter<unsigned int>("layersNeeded")),
0024 digitizeBadChambers_(p.getParameter<bool>("digitizeBadChambers")) {
0025 if (p.getParameter<bool>("doNeutrons")) {
0026 theNeutronReader = new CSCNeutronReader(p.getParameter<edm::ParameterSet>("neutrons"));
0027 }
0028 }
0029
0030 CSCDigitizer::~CSCDigitizer() {
0031 delete theNeutronReader;
0032 delete theStripElectronicsSim;
0033 delete theWireElectronicsSim;
0034 delete theStripHitSim;
0035 delete theWireHitSim;
0036 delete theDriftSim;
0037 }
0038
0039 void CSCDigitizer::doAction(MixCollection<PSimHit> &simHits,
0040 CSCWireDigiCollection &wireDigis,
0041 CSCStripDigiCollection &stripDigis,
0042 CSCComparatorDigiCollection &comparators,
0043 DigiSimLinks &wireDigiSimLinks,
0044 DigiSimLinks &stripDigiSimLinks,
0045 CLHEP::HepRandomEngine *engine) {
0046
0047 std::map<int, edm::PSimHitContainer> hitMap;
0048 for (MixCollection<PSimHit>::MixItr hitItr = simHits.begin(); hitItr != simHits.end(); ++hitItr) {
0049 hitMap[hitItr->detUnitId()].push_back(*hitItr);
0050 }
0051
0052
0053 std::map<int, std::set<int>> layersInChamberHit;
0054 for (std::map<int, edm::PSimHitContainer>::const_iterator hitMapItr = hitMap.begin(); hitMapItr != hitMap.end();
0055 ++hitMapItr) {
0056 CSCDetId cscDetId(hitMapItr->first);
0057 int chamberId = cscDetId.chamberId();
0058 layersInChamberHit[chamberId].insert(cscDetId.layer());
0059 }
0060
0061
0062 if (theNeutronReader != nullptr) {
0063 theNeutronReader->addHits(hitMap, engine);
0064 }
0065
0066
0067 for (std::map<int, edm::PSimHitContainer>::const_iterator hitMapItr = hitMap.begin(); hitMapItr != hitMap.end();
0068 ++hitMapItr) {
0069 CSCDetId detId = CSCDetId(hitMapItr->first);
0070 int chamberId = detId.chamberId();
0071 int endc = detId.endcap();
0072 int stat = detId.station();
0073 int ring = detId.ring();
0074 int cham = detId.chamber();
0075
0076 unsigned int nLayersInChamberHitForWireDigis = 0;
0077 if (stat == 1 && ring == 1) {
0078 std::set<int> layersInME1a = layersInChamberHit[CSCDetId(endc, stat, 4, cham, 0)];
0079 std::set<int> layersInME11 = layersInChamberHit[chamberId];
0080 layersInME11.insert(layersInME1a.begin(), layersInME1a.end());
0081 nLayersInChamberHitForWireDigis = layersInME11.size();
0082 } else if (stat == 1 && ring == 4) {
0083 std::set<int> layersInME1b = layersInChamberHit[CSCDetId(endc, stat, 1, cham, 0)];
0084 std::set<int> layersInME11 = layersInChamberHit[chamberId];
0085 layersInME11.insert(layersInME1b.begin(), layersInME1b.end());
0086 nLayersInChamberHitForWireDigis = layersInME11.size();
0087 } else
0088 nLayersInChamberHitForWireDigis = layersInChamberHit[chamberId].size();
0089
0090 unsigned int nLayersInChamberHitForStripDigis = layersInChamberHit[chamberId].size();
0091
0092 if (nLayersInChamberHitForWireDigis < theLayersNeeded && nLayersInChamberHitForStripDigis < theLayersNeeded)
0093 continue;
0094
0095 if (!digitizeBadChambers_ && theConditions->isInBadChamber(CSCDetId(hitMapItr->first)))
0096 continue;
0097
0098 const CSCLayer *layer = findLayer(hitMapItr->first);
0099 const edm::PSimHitContainer &layerSimHits = hitMapItr->second;
0100
0101 std::vector<CSCDetectorHit> newWireHits, newStripHits;
0102
0103
0104
0105 edm::LogVerbatim("CSCDigitizer") << "[CSCDigitizer] found " << layerSimHits.size() << " hit(s) in layer"
0106 << layer->id();
0107
0108
0109
0110
0111
0112 { newWireHits.swap(theWireHitSim->simulate(layer, layerSimHits, engine)); }
0113 if (!newWireHits.empty()) {
0114 newStripHits.swap(theStripHitSim->simulate(layer, newWireHits));
0115 }
0116
0117
0118 if (nLayersInChamberHitForWireDigis >= theLayersNeeded) {
0119 theWireElectronicsSim->simulate(layer, newWireHits, engine);
0120 theWireElectronicsSim->fillDigis(wireDigis, engine);
0121 wireDigiSimLinks.insert(theWireElectronicsSim->digiSimLinks());
0122 }
0123 if (nLayersInChamberHitForStripDigis >= theLayersNeeded) {
0124 theStripElectronicsSim->simulate(layer, newStripHits, engine);
0125 theStripElectronicsSim->fillDigis(stripDigis, comparators, engine);
0126 stripDigiSimLinks.insert(theStripElectronicsSim->digiSimLinks());
0127 }
0128 }
0129
0130
0131 std::list<int> missingLayers = layersMissing(stripDigis);
0132 for (std::list<int>::const_iterator missingLayerItr = missingLayers.begin(); missingLayerItr != missingLayers.end();
0133 ++missingLayerItr) {
0134 const CSCLayer *layer = findLayer(*missingLayerItr);
0135 theStripElectronicsSim->fillMissingLayer(layer, comparators, stripDigis, engine);
0136 }
0137 }
0138
0139 std::list<int> CSCDigitizer::layersMissing(const CSCStripDigiCollection &stripDigis) const {
0140 std::list<int> result;
0141
0142 std::map<int, std::list<int>> layersInChamberWithDigi;
0143 for (CSCStripDigiCollection::DigiRangeIterator j = stripDigis.begin(); j != stripDigis.end(); j++) {
0144 CSCDetId layerId((*j).first);
0145
0146 if ((*j).second.first != (*j).second.second) {
0147 int chamberId = layerId.chamberId();
0148 layersInChamberWithDigi[chamberId].push_back(layerId.layer());
0149 }
0150 }
0151
0152 std::list<int> oneThruSix;
0153 for (int i = 1; i <= 6; ++i)
0154 oneThruSix.push_back(i);
0155
0156 for (std::map<int, std::list<int>>::iterator layersInChamberWithDigiItr = layersInChamberWithDigi.begin();
0157 layersInChamberWithDigiItr != layersInChamberWithDigi.end();
0158 ++layersInChamberWithDigiItr) {
0159 std::list<int> &layersHit = layersInChamberWithDigiItr->second;
0160 if (layersHit.size() < 6 && layersHit.size() >= theLayersNeeded) {
0161 layersHit.sort();
0162 std::list<int> missingLayers(6);
0163 std::list<int>::iterator lastLayerMissing = set_difference(
0164 oneThruSix.begin(), oneThruSix.end(), layersHit.begin(), layersHit.end(), missingLayers.begin());
0165 int chamberId = layersInChamberWithDigiItr->first;
0166 for (std::list<int>::iterator layerMissingItr = missingLayers.begin(); layerMissingItr != lastLayerMissing;
0167 ++layerMissingItr) {
0168
0169 result.push_back(chamberId + *layerMissingItr);
0170 }
0171 }
0172 }
0173 return result;
0174 }
0175
0176 void CSCDigitizer::setMagneticField(const MagneticField *field) { theDriftSim->setMagneticField(field); }
0177
0178 void CSCDigitizer::setStripConditions(CSCStripConditions *cond) {
0179 theConditions = cond;
0180 theStripElectronicsSim->setStripConditions(cond);
0181 }
0182
0183 void CSCDigitizer::setParticleDataTable(const ParticleDataTable *pdt) { theWireHitSim->setParticleDataTable(pdt); }
0184
0185 const CSCLayer *CSCDigitizer::findLayer(int detId) const {
0186 assert(theCSCGeometry != nullptr);
0187 const GeomDetUnit *detUnit = theCSCGeometry->idToDetUnit(CSCDetId(detId));
0188 if (detUnit == nullptr) {
0189 throw cms::Exception("CSCDigiProducer") << "Invalid DetUnit: " << CSCDetId(detId)
0190 << "\nPerhaps your signal or pileup dataset are not compatible with "
0191 "the current release?";
0192 }
0193 return dynamic_cast<const CSCLayer *>(detUnit);
0194 }