DTCalibration/src/DTSegmentSelector.cc

0001 #include "CalibMuon/DTCalibration/interface/DTSegmentSelector.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Framework/interface/ESHandle.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
0010 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
0011 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0012 #include "DataFormats/DTRecHit/interface/DTRecHit1D.h"
0013 #include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
0014 #include "CondFormats/DTObjects/interface/DTStatusFlag.h"
0015
0016 DTSegmentSelector::DTSegmentSelector(edm::ParameterSet const& pset, edm::ConsumesCollector& iC)
0017     : muonTags_(pset.getParameter<edm::InputTag>("Muons")),
0018       checkNoisyChannels_(pset.getParameter<bool>("checkNoisyChannels")),
0019       minHitsPhi_(pset.getParameter<int>("minHitsPhi")),
0020       minHitsZ_(pset.getParameter<int>("minHitsZ")),
0021       maxChi2_(pset.getParameter<double>("maxChi2")),
0022       maxAnglePhi_(pset.getParameter<double>("maxAnglePhi")),
0023       maxAngleZ_(pset.getParameter<double>("maxAngleZ")) {
0024   muonToken_ = iC.consumes<reco::MuonCollection>(muonTags_);
0025   theStatusMapToken_ = iC.esConsumes();
0026 }
0027
0028 bool DTSegmentSelector::operator()(DTRecSegment4D const& segment,
0029                                    edm::Event const& event,
0030                                    edm::EventSetup const& setup) {
0031   bool result = true;
0032
0033   /* get the muon collection if one is specified
0034     (not specifying a muon collection switches off muon matching */
0035   if (!muonTags_.label().empty()) {
0036     edm::Handle<reco::MuonCollection> muonH;
0037     event.getByToken(muonToken_, muonH);
0038     const std::vector<reco::Muon>* muons = muonH.product();
0039     //std::cout << " Muon collection size: " << muons->size() << std::endl;
0040     if (muons->empty())
0041       return false;
0042
0043     DTChamberId segId(segment.rawId());
0044
0045     bool matched = false;
0046     for (auto& imuon : *muons)
0047       for (const auto& ch : imuon.matches()) {
0048         DetId chId(ch.id.rawId());
0049         if (chId != segId)
0050           continue;
0051         if (imuon.pt() < 15 || !imuon.isGlobalMuon())
0052           continue;
0053
0054         int nsegs = ch.segmentMatches.size();
0055         if (!nsegs)
0056           continue;
0057
0058         LocalPoint posHit = segment.localPosition();
0059         float dx = (posHit.x() ? posHit.x() - ch.x : 0);
0060         float dy = (posHit.y() ? posHit.y() - ch.y : 0);
0061         float dr = sqrt(dx * dx + dy * dy);
0062         if (dr < 5)
0063           matched = true;
0064       }
0065
0066     if (!matched)
0067       result = false;
0068   }
0069
0070   edm::ESHandle<DTStatusFlag> statusMap;
0071   if (checkNoisyChannels_)
0072     statusMap = setup.getHandle(theStatusMapToken_);
0073
0074   // Get the Phi 2D segment
0075   int nPhiHits = -1;
0076   bool segmentNoisyPhi = false;
0077   if (segment.hasPhi()) {
0078     const DTChamberRecSegment2D* phiSeg = segment.phiSegment();  // phiSeg lives in the chamber RF
0079     std::vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
0080     nPhiHits = phiRecHits.size();
0081     if (checkNoisyChannels_)
0082       segmentNoisyPhi = checkNoisySegment(statusMap, phiRecHits);
0083     //} else result = false;
0084   }
0085   // Get the Theta 2D segment
0086   int nZHits = -1;
0087   bool segmentNoisyZ = false;
0088   if (segment.hasZed()) {
0089     const DTSLRecSegment2D* zSeg = segment.zSegment();  // zSeg lives in the SL RF
0090     std::vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
0091     nZHits = zRecHits.size();
0092     if (checkNoisyChannels_)
0093       segmentNoisyZ = checkNoisySegment(statusMap, zRecHits);
0094   }
0095
0096   // Segment selection
0097   // Discard segment if it has a noisy cell
0098   if (segmentNoisyPhi || segmentNoisyZ)
0099     result = false;
0100
0101   // 2D-segment number of hits
0102   if (segment.hasPhi() && nPhiHits < minHitsPhi_)
0103     result = false;
0104
0105   if (segment.hasZed() && nZHits < minHitsZ_)
0106     result = false;
0107
0108   // Segment chi2
0109   double chiSquare = segment.chi2() / segment.degreesOfFreedom();
0110   if (chiSquare > maxChi2_)
0111     result = false;
0112
0113   // Segment angle
0114   LocalVector segment4DLocalDir = segment.localDirection();
0115   double angleZ = fabs(atan(segment4DLocalDir.y() / segment4DLocalDir.z()) * 180. / Geom::pi());
0116   if (angleZ > maxAngleZ_)
0117     result = false;
0118
0119   double anglePhi = fabs(atan(segment4DLocalDir.x() / segment4DLocalDir.z()) * 180. / Geom::pi());
0120   if (anglePhi > maxAnglePhi_)
0121     result = false;
0122
0123   return result;
0124 }
0125
0126 bool DTSegmentSelector::checkNoisySegment(edm::ESHandle<DTStatusFlag> const& statusMap,
0127                                           std::vector<DTRecHit1D> const& dtHits) {
0128   bool segmentNoisy = false;
0129
0130   std::vector<DTRecHit1D>::const_iterator dtHit = dtHits.begin();
0131   std::vector<DTRecHit1D>::const_iterator dtHits_end = dtHits.end();
0132   for (; dtHit != dtHits_end; ++dtHit) {
0133     DTWireId wireId = dtHit->wireId();
0134     // Check for noisy channels to skip them
0135     bool isNoisy = false, isFEMasked = false, isTDCMasked = false, isTrigMask = false, isDead = false, isNohv = false;
0136     statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
0137     if (isNoisy) {
0138       LogTrace("Calibration") << "Wire: " << wireId << " is noisy, skipping!";
0139       segmentNoisy = true;
0140       break;
0141     }
0142   }
0143   return segmentNoisy;
0144 }