CSCDigitizer/src/CSCDigitizer.cc

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   // arrange the hits by layer
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   // count how many layers on each chamber are hit
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   // add neutron background, if needed
0062   if (theNeutronReader != nullptr) {
0063     theNeutronReader->addHits(hitMap, engine);
0064   }
0065
0066   // now loop over layers and run the simulation for each one
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) {  // ME1b
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) {  // ME1a
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     // skip bad chambers
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     //    LogTrace("CSCDigitizer") << "CSCDigitizer: found " <<
0104     //    layerSimHits.size() <<" hit(s) in layer"
0105     edm::LogVerbatim("CSCDigitizer") << "[CSCDigitizer] found " << layerSimHits.size() << " hit(s) in layer"
0106                                      << layer->id();
0107     //       << " E" << layer->id().endcap() << " S" << layer->id().station() <<
0108     //       " R" << layer->id().ring()
0109     //       << " C" << layer->id().chamber() << " L" << layer->id().layer();
0110
0111     // turn the edm::PSimHits into WireHits, using the WireHitSim
0112     { newWireHits.swap(theWireHitSim->simulate(layer, layerSimHits, engine)); }
0113     if (!newWireHits.empty()) {
0114       newStripHits.swap(theStripHitSim->simulate(layer, newWireHits));
0115     }
0116
0117     // turn the hits into wire digis, using the electronicsSim
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   // fill in the layers were missing from this chamber
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     // make sure the vector of digis isn't empty
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         // got from layer 1-6 to layer ID
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;                              // cache here
0180   theStripElectronicsSim->setStripConditions(cond);  // percolate downwards
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 }