MCTruth/src/TrackerMuonHitExtractor.cc

0001 //
0002 // modified & integrated by Giovanni Abbiendi
0003 // from code by Arun Luthra:
0004 // UserCode/luthra/MuonTrackSelector/src/MuonTrackSelector.cc
0005 //
0006 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0007 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0008 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0009 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0010 #include "DataFormats/TrackReco/interface/Track.h"
0011 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "SimMuon/MCTruth/interface/TrackerMuonHitExtractor.h"
0014 #include <sstream>
0015
0016 TrackerMuonHitExtractor::TrackerMuonHitExtractor(const edm::ParameterSet &parset, edm::ConsumesCollector &&ic)
0017     : inputDTRecSegment4DToken_(
0018           ic.consumes<DTRecSegment4DCollection>(parset.getParameter<edm::InputTag>("inputDTRecSegment4DCollection"))),
0019       inputCSCSegmentToken_(
0020           ic.consumes<CSCSegmentCollection>(parset.getParameter<edm::InputTag>("inputCSCSegmentCollection"))) {}
0021
0022 void TrackerMuonHitExtractor::init(const edm::Event &iEvent) {
0023   const edm::Handle<DTRecSegment4DCollection> &dtSegmentCollectionH = iEvent.getHandle(inputDTRecSegment4DToken_);
0024   const edm::Handle<CSCSegmentCollection> &cscSegmentCollectionH = iEvent.getHandle(inputCSCSegmentToken_);
0025
0026   edm::LogVerbatim("TrackerMuonHitExtractor") << "\nThere are " << dtSegmentCollectionH->size() << " DT segments.";
0027   unsigned int index_dt_segment = 0;
0028   for (DTRecSegment4DCollection::const_iterator segment = dtSegmentCollectionH->begin();
0029        segment != dtSegmentCollectionH->end();
0030        ++segment, index_dt_segment++) {
0031     LocalPoint segmentLocalPosition = segment->localPosition();
0032     LocalVector segmentLocalDirection = segment->localDirection();
0033     LocalError segmentLocalPositionError = segment->localPositionError();
0034     LocalError segmentLocalDirectionError = segment->localDirectionError();
0035     DetId geoid = segment->geographicalId();
0036     DTChamberId dtdetid = DTChamberId(geoid);
0037     int wheel = dtdetid.wheel();
0038     int station = dtdetid.station();
0039     int sector = dtdetid.sector();
0040
0041     float segmentX = segmentLocalPosition.x();
0042     float segmentY = segmentLocalPosition.y();
0043     float segmentdXdZ = segmentLocalDirection.x() / segmentLocalDirection.z();
0044     float segmentdYdZ = segmentLocalDirection.y() / segmentLocalDirection.z();
0045     float segmentXerr = sqrt(segmentLocalPositionError.xx());
0046     float segmentYerr = sqrt(segmentLocalPositionError.yy());
0047     float segmentdXdZerr = sqrt(segmentLocalDirectionError.xx());
0048     float segmentdYdZerr = sqrt(segmentLocalDirectionError.yy());
0049
0050     edm::LogVerbatim("TrackerMuonHitExtractor")
0051         << "\nDT segment index :" << index_dt_segment << "\nchamber Wh:" << wheel << ",St:" << station
0052         << ",Se:" << sector << "\nLocal Position (X,Y)=(" << segmentX << "," << segmentY << ") +/- (" << segmentXerr
0053         << "," << segmentYerr << "), "
0054         << "Local Direction (dXdZ,dYdZ)=(" << segmentdXdZ << "," << segmentdYdZ << ") +/- (" << segmentdXdZerr << ","
0055         << segmentdYdZerr << ")";
0056   }
0057
0058   edm::LogVerbatim("TrackerMuonHitExtractor") << "\nThere are " << cscSegmentCollectionH->size() << " CSC segments.";
0059   unsigned int index_csc_segment = 0;
0060   for (CSCSegmentCollection::const_iterator segment = cscSegmentCollectionH->begin();
0061        segment != cscSegmentCollectionH->end();
0062        ++segment, index_csc_segment++) {
0063     LocalPoint segmentLocalPosition = segment->localPosition();
0064     LocalVector segmentLocalDirection = segment->localDirection();
0065     LocalError segmentLocalPositionError = segment->localPositionError();
0066     LocalError segmentLocalDirectionError = segment->localDirectionError();
0067
0068     DetId geoid = segment->geographicalId();
0069     CSCDetId cscdetid = CSCDetId(geoid);
0070     int endcap = cscdetid.endcap();
0071     int station = cscdetid.station();
0072     int ring = cscdetid.ring();
0073     int chamber = cscdetid.chamber();
0074
0075     float segmentX = segmentLocalPosition.x();
0076     float segmentY = segmentLocalPosition.y();
0077     float segmentdXdZ = segmentLocalDirection.x() / segmentLocalDirection.z();
0078     float segmentdYdZ = segmentLocalDirection.y() / segmentLocalDirection.z();
0079     float segmentXerr = sqrt(segmentLocalPositionError.xx());
0080     float segmentYerr = sqrt(segmentLocalPositionError.yy());
0081     float segmentdXdZerr = sqrt(segmentLocalDirectionError.xx());
0082     float segmentdYdZerr = sqrt(segmentLocalDirectionError.yy());
0083
0084     edm::LogVerbatim("TrackerMuonHitExtractor")
0085         << "\nCSC segment index :" << index_csc_segment << "\nchamber Endcap:" << endcap << ",St:" << station
0086         << ",Ri:" << ring << ",Ch:" << chamber << "\nLocal Position (X,Y)=(" << segmentX << "," << segmentY << ") +/- ("
0087         << segmentXerr << "," << segmentYerr << "), "
0088         << "Local Direction (dXdZ,dYdZ)=(" << segmentdXdZ << "," << segmentdYdZ << ") +/- (" << segmentdXdZerr << ","
0089         << segmentdYdZerr << ")";
0090   }
0091 }
0092 std::vector<const TrackingRecHit *> TrackerMuonHitExtractor::getMuonHits(const reco::Muon &mu) const {
0093   std::vector<const TrackingRecHit *> ret;
0094
0095   int wheel, station, sector;
0096   int endcap, /*station, */ ring, chamber;
0097
0098   edm::LogVerbatim("TrackerMuonHitExtractor")
0099       << "Number of chambers: " << mu.matches().size()
0100       << ", arbitrated: " << mu.numberOfMatches(reco::Muon::SegmentAndTrackArbitration)
0101       << ", tot (no arbitration): " << mu.numberOfMatches(reco::Muon::NoArbitration);
0102   unsigned int index_chamber = 0;
0103   int n_segments_noArb = 0;
0104
0105   for (std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = mu.matches().begin();
0106        chamberMatch != mu.matches().end();
0107        ++chamberMatch, index_chamber++) {
0108     std::stringstream chamberStr;
0109     chamberStr << "\nchamber index: " << index_chamber;
0110
0111     int subdet = chamberMatch->detector();
0112     DetId did = chamberMatch->id;
0113
0114     if (subdet == MuonSubdetId::DT) {
0115       DTChamberId dtdetid = DTChamberId(did);
0116       wheel = dtdetid.wheel();
0117       station = dtdetid.station();
0118       sector = dtdetid.sector();
0119       chamberStr << ", DT chamber Wh:" << wheel << ",St:" << station << ",Se:" << sector;
0120     } else if (subdet == MuonSubdetId::CSC) {
0121       CSCDetId cscdetid = CSCDetId(did);
0122       endcap = cscdetid.endcap();
0123       station = cscdetid.station();
0124       ring = cscdetid.ring();
0125       chamber = cscdetid.chamber();
0126       chamberStr << ", CSC chamber End:" << endcap << ",St:" << station << ",Ri:" << ring << ",Ch:" << chamber;
0127     }
0128
0129     chamberStr << ", Number of segments: " << chamberMatch->segmentMatches.size();
0130     edm::LogVerbatim("TrackerMuonHitExtractor") << chamberStr.str();
0131     n_segments_noArb = n_segments_noArb + chamberMatch->segmentMatches.size();
0132
0133     unsigned int index_segment = 0;
0134
0135     for (std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();
0136          segmentMatch != chamberMatch->segmentMatches.end();
0137          ++segmentMatch, index_segment++) {
0138       float segmentX = segmentMatch->x;
0139       float segmentY = segmentMatch->y;
0140       float segmentdXdZ = segmentMatch->dXdZ;
0141       float segmentdYdZ = segmentMatch->dYdZ;
0142       float segmentXerr = segmentMatch->xErr;
0143       float segmentYerr = segmentMatch->yErr;
0144       float segmentdXdZerr = segmentMatch->dXdZErr;
0145       float segmentdYdZerr = segmentMatch->dYdZErr;
0146
0147       CSCSegmentRef segmentCSC = segmentMatch->cscSegmentRef;
0148       DTRecSegment4DRef segmentDT = segmentMatch->dtSegmentRef;
0149
0150       bool segment_arbitrated_Ok = (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) &&
0151                                     segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR));
0152
0153       std::string ARBITRATED(" ***Arbitrated Off*** ");
0154       if (segment_arbitrated_Ok)
0155         ARBITRATED = " ***ARBITRATED OK*** ";
0156
0157       if (subdet == MuonSubdetId::DT) {
0158         edm::LogVerbatim("TrackerMuonHitExtractor")
0159             << "\n\t segment index: " << index_segment << ARBITRATED << "\n\t  Local Position (X,Y)=(" << segmentX
0160             << "," << segmentY << ") +/- (" << segmentXerr << "," << segmentYerr << "), "
0161             << "\n\t  Local Direction (dXdZ,dYdZ)=(" << segmentdXdZ << "," << segmentdYdZ << ") +/- (" << segmentdXdZerr
0162             << "," << segmentdYdZerr << ")";
0163
0164         // with the following line the DT segments failing standard arbitration
0165         // are skipped
0166         if (!segment_arbitrated_Ok)
0167           continue;
0168
0169         if (segmentDT.get() != nullptr) {
0170           const DTRecSegment4D *segment = segmentDT.get();
0171
0172           edm::LogVerbatim("TrackerMuonHitExtractor")
0173               << "\t ===> MATCHING with DT segment with index = " << segmentDT.key();
0174
0175           if (segment->hasPhi()) {
0176             const DTChamberRecSegment2D *phiSeg = segment->phiSegment();
0177             std::vector<const TrackingRecHit *> phiHits = phiSeg->recHits();
0178             for (std::vector<const TrackingRecHit *>::const_iterator ihit = phiHits.begin(); ihit != phiHits.end();
0179                  ++ihit) {
0180               ret.push_back(*ihit);
0181             }
0182           }
0183
0184           if (segment->hasZed()) {
0185             const DTSLRecSegment2D *zSeg = (*segment).zSegment();
0186             std::vector<const TrackingRecHit *> zedHits = zSeg->recHits();
0187             for (std::vector<const TrackingRecHit *>::const_iterator ihit = zedHits.begin(); ihit != zedHits.end();
0188                  ++ihit) {
0189               ret.push_back(*ihit);
0190             }
0191           }
0192         } else
0193           edm::LogWarning("TrackerMuonHitExtractor") << "\n***WARNING: UNMATCHED DT segment ! \n";
0194       }  // if (subdet == MuonSubdetId::DT)
0195
0196       else if (subdet == MuonSubdetId::CSC) {
0197         edm::LogVerbatim("TrackerMuonHitExtractor")
0198             << "\n\t segment index: " << index_segment << ARBITRATED << "\n\t  Local Position (X,Y)=(" << segmentX
0199             << "," << segmentY << ") +/- (" << segmentXerr << "," << segmentYerr << "), "
0200             << "\n\t  Local Direction (dXdZ,dYdZ)=(" << segmentdXdZ << "," << segmentdYdZ << ") +/- (" << segmentdXdZerr
0201             << "," << segmentdYdZerr << ")";
0202
0203         // with the following line the CSC segments failing standard arbitration
0204         // are skipped
0205         if (!segment_arbitrated_Ok)
0206           continue;
0207
0208         if (segmentCSC.get() != nullptr) {
0209           const CSCSegment *segment = segmentCSC.get();
0210
0211           edm::LogVerbatim("TrackerMuonHitExtractor")
0212               << "\t ===> MATCHING with CSC segment with index = " << segmentCSC.key();
0213
0214           std::vector<const TrackingRecHit *> hits = segment->recHits();
0215           for (std::vector<const TrackingRecHit *>::const_iterator ihit = hits.begin(); ihit != hits.end(); ++ihit) {
0216             ret.push_back(*ihit);
0217           }
0218         } else
0219           edm::LogWarning("TrackerMuonHitExtractor") << "\n***WARNING: UNMATCHED CSC segment ! \n";
0220       }  // else if (subdet == MuonSubdetId::CSC)
0221
0222     }  // loop on vector<MuonSegmentMatch>
0223   }  // loop on vector<MuonChamberMatch>
0224
0225   edm::LogVerbatim("TrackerMuonHitExtractor") << "\n N. matched Segments before arbitration = " << n_segments_noArb;
0226
0227   return ret;
0228 }