MCEmbeddingTools/plugins/MuonDetCleaner.cc

0001 #include "TauAnalysis/MCEmbeddingTools/plugins/MuonDetCleaner.h"
0002
0003 #include "DataFormats/DTRecHit/interface/DTRecHit1DPair.h"
0004 #include "DataFormats/DTRecHit/interface/DTSLRecCluster.h"
0005 #include "DataFormats/MuonDetId/interface/DTLayerId.h"
0006
0007 #include "DataFormats/CSCRecHit/interface/CSCRecHit2D.h"
0008 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0009 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0010 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
0011
0012 typedef MuonDetCleaner<CSCDetId, CSCRecHit2D> CSCRecHitColCleaner;
0013 typedef MuonDetCleaner<CSCDetId, CSCSegment> CSCSegmentColCleaner;
0014 typedef MuonDetCleaner<DTChamberId, DTRecSegment4D> DTRecSegment4DColCleaner;
0015 typedef MuonDetCleaner<DTLayerId, DTRecHit1DPair> DTRecHitColCleaner;
0016 typedef MuonDetCleaner<RPCDetId, RPCRecHit> RPCRecHitColCleaner;
0017
0018 template <typename T1, typename T2>
0019 MuonDetCleaner<T1, T2>::MuonDetCleaner(const edm::ParameterSet &iConfig)
0020     : mu_input_(consumes<edm::View<pat::Muon>>(iConfig.getParameter<edm::InputTag>("MuonCollection"))),
0021       m_dtDigisToken(consumes<DTDigiCollection>(iConfig.getParameter<edm::InputTag>("dtDigiCollectionLabel"))),
0022       m_cscDigisToken(consumes<CSCStripDigiCollection>(iConfig.getParameter<edm::InputTag>("cscDigiCollectionLabel"))),
0023       m_dtGeometryToken(esConsumes<edm::Transition::BeginRun>()),
0024       m_cscGeometryToken(esConsumes<edm::Transition::BeginRun>()),
0025       propagatorToken_(esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0026       m_digiMaxDistanceX(iConfig.getParameter<double>("digiMaxDistanceX")) {
0027   std::vector<edm::InputTag> inCollections = iConfig.getParameter<std::vector<edm::InputTag>>("oldCollection");
0028   for (const auto &inCollection : inCollections) {
0029     inputs_[inCollection.instance()] = consumes<RecHitCollection>(inCollection);
0030     produces<RecHitCollection>(inCollection.instance());
0031   }
0032   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0033   edm::ConsumesCollector iC = consumesCollector();
0034   parameters_.loadParameters(parameters, iC);
0035 }
0036
0037 template <typename T1, typename T2>
0038 MuonDetCleaner<T1, T2>::~MuonDetCleaner() {
0039   // nothing to be done yet...
0040 }
0041
0042 template <typename T1, typename T2>
0043 void MuonDetCleaner<T1, T2>::beginRun(const edm::Run &iRun, const edm::EventSetup &iSetup) {
0044   // get the geometries from the event setup everytime the run (of the event which is processed right now to the previous event) changes
0045   m_dtGeometry = iSetup.getHandle(m_dtGeometryToken);
0046   m_cscGeometry = iSetup.getHandle(m_cscGeometryToken);
0047 }
0048
0049 template <typename T1, typename T2>
0050 void MuonDetCleaner<T1, T2>::produce(edm::Event &iEvent, edm::EventSetup const &iSetup) {
0051   std::map<T1, std::vector<T2>> recHits_output;  // This data format is easyer to handle
0052   std::vector<uint32_t> vetoHits;
0053
0054   // First fill the veto RecHits colletion with the Hits from the input muons
0055   edm::Handle<edm::View<pat::Muon>> muonHandle;
0056   iEvent.getByToken(mu_input_, muonHandle);
0057   edm::View<pat::Muon> muons = *muonHandle;
0058   for (edm::View<pat::Muon>::const_iterator iMuon = muons.begin(); iMuon != muons.end(); ++iMuon) {
0059     const reco::Track *track = nullptr;
0060     if (iMuon->isGlobalMuon())
0061       track = iMuon->outerTrack().get();
0062     else if (iMuon->isStandAloneMuon())
0063       track = iMuon->outerTrack().get();
0064     else if (iMuon->isRPCMuon())
0065       track = iMuon->innerTrack().get();  // To add, try to access the rpc track
0066     else if (iMuon->isTrackerMuon())
0067       track = iMuon->innerTrack().get();
0068     else {
0069       edm::LogError("TauEmbedding") << "The imput muon: " << (*iMuon)
0070                                     << " must be either global or does or be tracker muon";
0071       assert(0);
0072     }
0073
0074     for (trackingRecHit_iterator hitIt = track->recHitsBegin(); hitIt != track->recHitsEnd(); ++hitIt) {
0075       const TrackingRecHit &murechit = **hitIt;  // Base class for all rechits
0076       if (!(murechit).isValid())
0077         continue;
0078       if (!checkrecHit(murechit))
0079         continue;                         // Check if the hit belongs to a specifc detector section
0080       fillVetoHits(murechit, &vetoHits);  // Go back to the very basic rechits
0081     }
0082
0083     sort(vetoHits.begin(), vetoHits.end());
0084     vetoHits.erase(unique(vetoHits.begin(), vetoHits.end()), vetoHits.end());
0085     iEvent.getByToken(m_dtDigisToken, m_dtDigis);
0086     iEvent.getByToken(m_cscDigisToken, m_cscDigis);
0087     edm::ESHandle<Propagator> propagator;
0088     trackAssociator_.setPropagator(&iSetup.getData(propagatorToken_));
0089     TrackDetMatchInfo info =
0090         trackAssociator_.associate(iEvent, iSetup, *track, parameters_, TrackDetectorAssociator::Any);
0091
0092     // inspired from Muon Identification algorithm: https://github.com/cms-sw/cmssw/blob/3b943c0dbbdf4494cd66064a5a147301f38af295/RecoMuon/MuonIdentification/plugins/MuonIdProducer.cc#L911
0093     // and the MuonShowerDigiFiller: https://github.com/cms-sw/cmssw/blob/29909891e150c9781f4ade2a6f7b5beb0bd67a6e/RecoMuon/MuonIdentification/interface/MuonShowerDigiFiller.h & https://github.com/cms-sw/cmssw/blob/29909891e150c9781f4ade2a6f7b5beb0bd67a6e/RecoMuon/MuonIdentification/src/MuonShowerDigiFiller.cc
0094     for (const auto &chamber : info.chambers) {
0095       if (chamber.id.subdetId() == MuonSubdetId::RPC)
0096         continue;
0097
0098       reco::MuonChamberMatch matchedChamber;
0099
0100       const auto &lErr = chamber.tState.localError();
0101       const auto &lPos = chamber.tState.localPosition();
0102       const auto &lDir = chamber.tState.localDirection();
0103       const auto &localError = lErr.positionError();
0104
0105       matchedChamber.x = lPos.x();
0106       matchedChamber.y = lPos.y();
0107       matchedChamber.xErr = sqrt(localError.xx());
0108       matchedChamber.yErr = sqrt(localError.yy());
0109       matchedChamber.dXdZ = lDir.z() != 0 ? lDir.x() / lDir.z() : 9999;
0110       matchedChamber.dYdZ = lDir.z() != 0 ? lDir.y() / lDir.z() : 9999;
0111
0112       // DANGEROUS - compiler cannot guaranty parameters ordering
0113       AlgebraicSymMatrix55 trajectoryCovMatrix = lErr.matrix();
0114       matchedChamber.dXdZErr = trajectoryCovMatrix(1, 1) > 0 ? sqrt(trajectoryCovMatrix(1, 1)) : 0;
0115       matchedChamber.dYdZErr = trajectoryCovMatrix(2, 2) > 0 ? sqrt(trajectoryCovMatrix(2, 2)) : 0;
0116
0117       matchedChamber.edgeX = chamber.localDistanceX;
0118       matchedChamber.edgeY = chamber.localDistanceY;
0119
0120       matchedChamber.id = chamber.id;
0121
0122       // DT chamber
0123       if (matchedChamber.detector() == MuonSubdetId::DT) {
0124         double xTrack = matchedChamber.x;
0125
0126         for (int sl = 1; sl <= DTChamberId::maxSuperLayerId; sl += 2) {
0127           for (int layer = 1; layer <= DTChamberId::maxLayerId; ++layer) {
0128             const DTLayerId layerId(DTChamberId(matchedChamber.id.rawId()), sl, layer);
0129             auto range = m_dtDigis->get(layerId);
0130
0131             for (auto digiIt = range.first; digiIt != range.second; ++digiIt) {
0132               const auto topo = m_dtGeometry->layer(layerId)->specificTopology();
0133               double xWire = topo.wirePosition((*digiIt).wire());
0134               double dX = std::abs(xWire - xTrack);
0135
0136               if (dX < m_digiMaxDistanceX) {
0137                 vetoHits.push_back(matchedChamber.id.rawId());
0138               }
0139             }
0140           }
0141         }
0142       }
0143
0144       else if (matchedChamber.detector() == MuonSubdetId::CSC) {
0145         double xTrack = matchedChamber.x;
0146         double yTrack = matchedChamber.y;
0147
0148         for (int iLayer = 1; iLayer <= CSCDetId::maxLayerId(); ++iLayer) {
0149           const CSCDetId chId(matchedChamber.id.rawId());
0150           const CSCDetId layerId(chId.endcap(), chId.station(), chId.ring(), chId.chamber(), iLayer);
0151           auto range = m_cscDigis->get(layerId);
0152
0153           for (auto digiIt = range.first; digiIt != range.second; ++digiIt) {
0154             std::vector<int> adcVals = digiIt->getADCCounts();
0155             bool hasFired = false;
0156             float pedestal = 0.5 * (float)(adcVals[0] + adcVals[1]);
0157             float threshold = 13.3;
0158             float diff = 0.;
0159             for (const auto &adcVal : adcVals) {
0160               diff = (float)adcVal - pedestal;
0161               if (diff > threshold) {
0162                 hasFired = true;
0163                 break;
0164               }
0165             }
0166
0167             if (!hasFired)
0168               continue;
0169
0170             const CSCLayerGeometry *layerGeom = m_cscGeometry->layer(layerId)->geometry();
0171             Float_t xStrip = layerGeom->xOfStrip(digiIt->getStrip(), yTrack);
0172             float dX = std::abs(xStrip - xTrack);
0173
0174             if (dX < m_digiMaxDistanceX) {
0175               vetoHits.push_back(matchedChamber.id.rawId());
0176             }
0177           }
0178         }
0179       }
0180     }
0181
0182     // End improved cleaning in the muon chambers
0183   }
0184
0185   sort(vetoHits.begin(), vetoHits.end());
0186   vetoHits.erase(unique(vetoHits.begin(), vetoHits.end()), vetoHits.end());
0187
0188   // Now this can also handle different instance
0189   for (auto input_ : inputs_) {
0190     // Second read in the RecHit Colltection which is to be replaced, without the vetoRecHits
0191     typedef edm::Handle<RecHitCollection> RecHitCollectionHandle;
0192     RecHitCollectionHandle RecHitinput;
0193     iEvent.getByToken(input_.second, RecHitinput);
0194     for (typename RecHitCollection::const_iterator recHit = RecHitinput->begin(); recHit != RecHitinput->end();
0195          ++recHit) {  // loop over the basic rec hit collection (DT CSC or RPC)
0196       if (find(vetoHits.begin(), vetoHits.end(), getRawDetId(*recHit)) != vetoHits.end())
0197         continue;  // If the hit is not in the
0198       T1 detId(getRawDetId(*recHit));
0199       recHits_output[detId].push_back(*recHit);
0200     }
0201
0202     // Last step savet the output in the CMSSW Data Format
0203     std::unique_ptr<RecHitCollection> output(new RecHitCollection());
0204     for (typename std::map<T1, std::vector<T2>>::const_iterator recHit = recHits_output.begin();
0205          recHit != recHits_output.end();
0206          ++recHit) {
0207       output->put(recHit->first, recHit->second.begin(), recHit->second.end());
0208     }
0209     output->post_insert();
0210     iEvent.put(std::move(output), input_.first);
0211   }
0212 }
0213
0214 template <typename T1, typename T2>
0215 void MuonDetCleaner<T1, T2>::fillVetoHits(const TrackingRecHit &rh, std::vector<uint32_t> *HitsList) {
0216   std::vector<const TrackingRecHit *> rh_components = rh.recHits();
0217   if (rh_components.empty()) {
0218     HitsList->push_back(rh.rawId());
0219   } else {
0220     for (std::vector<const TrackingRecHit *>::const_iterator rh_component = rh_components.begin();
0221          rh_component != rh_components.end();
0222          ++rh_component) {
0223       fillVetoHits(**rh_component, HitsList);
0224     }
0225   }
0226 }
0227
0228 //-------------------------------------------------------------------------------
0229 // define 'getRawDetId' functions used for different types of recHits
0230 //-------------------------------------------------------------------------------
0231
0232 template <typename T1, typename T2>
0233 uint32_t MuonDetCleaner<T1, T2>::getRawDetId(const T2 &recHit) {
0234   assert(0);  // CV: make sure general function never gets called;
0235               //     always use template specializations
0236 }
0237
0238 template <>
0239 uint32_t MuonDetCleaner<CSCDetId, CSCRecHit2D>::getRawDetId(const CSCRecHit2D &recHit) {
0240   return recHit.cscDetId().rawId();
0241 }
0242
0243 template <>
0244 uint32_t MuonDetCleaner<DTLayerId, DTRecHit1DPair>::getRawDetId(const DTRecHit1DPair &recHit) {
0245   return recHit.geographicalId().rawId();
0246 }
0247
0248 template <>
0249 uint32_t MuonDetCleaner<RPCDetId, RPCRecHit>::getRawDetId(const RPCRecHit &recHit) {
0250   return recHit.rpcId().rawId();
0251 }
0252
0253 //-------------------------------------------------------------------------------
0254 // find out what the kind of RecHit used by imput muons rechit
0255 //-------------------------------------------------------------------------------
0256
0257 template <typename T1, typename T2>
0258 bool MuonDetCleaner<T1, T2>::checkrecHit(const TrackingRecHit &recHit) {
0259   edm::LogError("TauEmbedding") << "!!!! Please add the checkrecHit for the individual class templates " assert(0);
0260 }
0261
0262 template <>
0263 bool MuonDetCleaner<CSCDetId, CSCRecHit2D>::checkrecHit(const TrackingRecHit &recHit) {
0264   const std::type_info &hit_type = typeid(recHit);
0265   if (hit_type == typeid(CSCSegment)) {
0266     return true;
0267   }  // This should be the default one (which are included in the global (outer) muon track)
0268   else if (hit_type == typeid(CSCRecHit2D)) {
0269     return true;
0270   }
0271   return false;
0272 }
0273
0274 template <>
0275 uint32_t MuonDetCleaner<CSCDetId, CSCSegment>::getRawDetId(const CSCSegment &recHit) {
0276   return recHit.cscDetId().rawId();
0277 }
0278
0279 template <>
0280 uint32_t MuonDetCleaner<DTChamberId, DTRecSegment4D>::getRawDetId(const DTRecSegment4D &recHit) {
0281   return recHit.geographicalId().rawId();
0282 }
0283
0284 template <>
0285 bool MuonDetCleaner<DTLayerId, DTRecHit1DPair>::checkrecHit(const TrackingRecHit &recHit) {
0286   const std::type_info &hit_type = typeid(recHit);
0287   if (hit_type == typeid(DTRecSegment4D)) {
0288     return true;
0289   }  // This should be the default one (which are included in the global (outer) muon track)
0290   else if (hit_type == typeid(DTRecHit1D)) {
0291     return true;
0292   } else if (hit_type == typeid(DTSLRecCluster)) {
0293     return true;
0294   } else if (hit_type == typeid(DTSLRecSegment2D)) {
0295     return true;
0296   }
0297   return false;
0298 }
0299
0300 template <>
0301 bool MuonDetCleaner<RPCDetId, RPCRecHit>::checkrecHit(const TrackingRecHit &recHit) {
0302   const std::type_info &hit_type = typeid(recHit);
0303   if (hit_type == typeid(RPCRecHit)) {
0304     return true;
0305   }  // This should be the default one (which are included in the global (outer) muon track)
0306   return false;
0307 }
0308
0309 template <>
0310 bool MuonDetCleaner<CSCDetId, CSCSegment>::checkrecHit(const TrackingRecHit &recHit) {
0311   const std::type_info &hit_type = typeid(recHit);
0312   if (hit_type == typeid(CSCSegment)) {
0313     return true;
0314   }  // This should be the default one (which are included in the global (outer) muon track)
0315   return false;
0316 }
0317
0318 template <>
0319 bool MuonDetCleaner<DTChamberId, DTRecSegment4D>::checkrecHit(const TrackingRecHit &recHit) {
0320   const std::type_info &hit_type = typeid(recHit);
0321   if (hit_type == typeid(DTRecSegment4D)) {
0322     return true;
0323   }  // This should be the default one (which are included in the global (outer) muon track)
0324   else if (hit_type == typeid(DTRecHit1D)) {
0325     return true;
0326   } else if (hit_type == typeid(DTSLRecCluster)) {
0327     return true;
0328   } else if (hit_type == typeid(DTSLRecSegment2D)) {
0329     return true;
0330   }
0331   return false;
0332 }
0333
0334 DEFINE_FWK_MODULE(CSCRecHitColCleaner);
0335 DEFINE_FWK_MODULE(DTRecHitColCleaner);
0336 DEFINE_FWK_MODULE(RPCRecHitColCleaner);
0337 DEFINE_FWK_MODULE(CSCSegmentColCleaner);
0338 DEFINE_FWK_MODULE(DTRecSegment4DColCleaner);