File indexing completed on 2023-03-17 11:20:35
0001 #include "RecoMuon/MuonIdentification/interface/MuonIdTruthInfo.h"
0002 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0003 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0004 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0005 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0008 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010
0011 void MuonIdTruthInfo::registerConsumes(edm::ConsumesCollector& iC) {
0012 iC.mayConsume<edm::SimTrackContainer>(edm::InputTag("g4SimHits", ""));
0013 iC.mayConsume<edm::PSimHitContainer>(edm::InputTag("g4SimHits", "MuonDTHits"));
0014 iC.mayConsume<edm::PSimHitContainer>(edm::InputTag("g4SimHits", "MuonCSCHits"));
0015 }
0016
0017 void MuonIdTruthInfo::truthMatchMuon(const edm::Event& iEvent,
0018 GlobalTrackingGeometry const& geometry,
0019 reco::Muon& aMuon) {
0020
0021
0022
0023
0024
0025 edm::Handle<edm::SimTrackContainer> simTracks;
0026 iEvent.getByLabel<edm::SimTrackContainer>("g4SimHits", "", simTracks);
0027 if (!simTracks.isValid()) {
0028 LogTrace("MuonIdentification") << "No tracks found";
0029 return;
0030 }
0031
0032 float bestMatchChi2 = 9999;
0033 unsigned int bestMatch = 0;
0034 const unsigned int offset = 0;
0035
0036 for (edm::SimTrackContainer::const_iterator simTrk = simTracks->begin(); simTrk != simTracks->end(); simTrk++) {
0037 float chi2 = matchChi2(*aMuon.track().get(), *simTrk);
0038 if (chi2 > bestMatchChi2)
0039 continue;
0040 bestMatchChi2 = chi2;
0041 bestMatch = simTrk->trackId();
0042 }
0043
0044 bestMatch -= offset;
0045
0046 std::vector<reco::MuonChamberMatch>& matches = aMuon.matches();
0047 int numberOfTruthMatchedChambers = 0;
0048
0049
0050 for (std::vector<reco::MuonChamberMatch>::iterator chamberMatch = matches.begin(); chamberMatch != matches.end();
0051 chamberMatch++) {
0052 if (chamberMatch->id.det() != DetId::Muon) {
0053 edm::LogWarning("MuonIdentification") << "Detector id of a muon chamber corresponds to not a muon detector";
0054 continue;
0055 }
0056 reco::MuonSegmentMatch bestSegmentMatch;
0057 double distance = 99999;
0058
0059 if (chamberMatch->id.subdetId() == MuonSubdetId::DT) {
0060 DTChamberId detId(chamberMatch->id.rawId());
0061
0062 edm::Handle<edm::PSimHitContainer> simHits;
0063 iEvent.getByLabel("g4SimHits", "MuonDTHits", simHits);
0064 if (simHits.isValid()) {
0065 for (edm::PSimHitContainer::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++)
0066 if (hit->trackId() == bestMatch)
0067 checkSimHitForBestMatch(bestSegmentMatch, distance, *hit, detId, geometry);
0068 } else
0069 LogTrace("MuonIdentification") << "No DT simulated hits are found";
0070 }
0071
0072 if (chamberMatch->id.subdetId() == MuonSubdetId::CSC) {
0073 CSCDetId detId(chamberMatch->id.rawId());
0074
0075 edm::Handle<edm::PSimHitContainer> simHits;
0076 iEvent.getByLabel("g4SimHits", "MuonCSCHits", simHits);
0077 if (simHits.isValid()) {
0078 for (edm::PSimHitContainer::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++)
0079 if (hit->trackId() == bestMatch)
0080 checkSimHitForBestMatch(bestSegmentMatch, distance, *hit, detId, geometry);
0081 } else
0082 LogTrace("MuonIdentification") << "No CSC simulated hits are found";
0083 }
0084 if (distance < 9999) {
0085 chamberMatch->truthMatches.push_back(bestSegmentMatch);
0086 numberOfTruthMatchedChambers++;
0087 LogTrace("MuonIdentification") << "Best truth matched hit:"
0088 << "\tDetId: " << chamberMatch->id.rawId() << "\n"
0089 << "\tprojection: ( " << bestSegmentMatch.x << ", " << bestSegmentMatch.y
0090 << " )\n";
0091 }
0092 }
0093 LogTrace("MuonIdentification") << "Truth matching summary:\n\tnumber of chambers: " << matches.size()
0094 << "\n\tnumber of truth matched chambers: " << numberOfTruthMatchedChambers << "\n";
0095 }
0096
0097 void MuonIdTruthInfo::checkSimHitForBestMatch(reco::MuonSegmentMatch& segmentMatch,
0098 double& distance,
0099 const PSimHit& hit,
0100 const DetId& chamberId,
0101 const GlobalTrackingGeometry& geometry) {
0102 printf("DONT FORGET TO CALL REGISTERCONSUMES()\n");
0103
0104
0105
0106
0107
0108
0109 const GeomDet* chamberGeometry = geometry.idToDet(chamberId);
0110 const GeomDet* simUnitGeometry = geometry.idToDet(DetId(hit.detUnitId()));
0111
0112 if (chamberGeometry && simUnitGeometry) {
0113 LocalPoint entryPoint = chamberGeometry->toLocal(simUnitGeometry->toGlobal(hit.entryPoint()));
0114 LocalPoint exitPoint = chamberGeometry->toLocal(simUnitGeometry->toGlobal(hit.exitPoint()));
0115 LocalVector direction = exitPoint - entryPoint;
0116 if (fabs(direction.z()) > 0.001) {
0117 LocalPoint projection = entryPoint - direction * (entryPoint.z() / direction.z());
0118 if (fabs(projection.z()) > 0.001)
0119 edm::LogWarning("MuonIdentification") << "z coordinate of the hit projection must be zero and it's not!\n";
0120
0121 double new_distance = 99999;
0122 if (entryPoint.z() * exitPoint.z() < -1)
0123 new_distance = 0;
0124 else {
0125 if (fabs(entryPoint.z()) < fabs(exitPoint.z()))
0126 new_distance = fabs(entryPoint.z());
0127 else
0128 new_distance = fabs(exitPoint.z());
0129 }
0130
0131 if (new_distance < distance) {
0132
0133 segmentMatch.x = projection.x();
0134 segmentMatch.y = projection.y();
0135 segmentMatch.xErr = 0;
0136 segmentMatch.yErr = 0;
0137 segmentMatch.dXdZ = direction.x() / direction.z();
0138 segmentMatch.dYdZ = direction.y() / direction.z();
0139 segmentMatch.dXdZErr = 0;
0140 segmentMatch.dYdZErr = 0;
0141 distance = new_distance;
0142 LogTrace("MuonIdentificationVerbose")
0143 << "Better truth matched segment found:\n"
0144 << "\tDetId: " << chamberId.rawId() << "\n"
0145 << "\tentry point: ( " << entryPoint.x() << ", " << entryPoint.y() << ", " << entryPoint.z() << " )\n"
0146 << "\texit point: ( " << exitPoint.x() << ", " << exitPoint.y() << ", " << exitPoint.z() << " )\n"
0147 << "\tprojection: ( " << projection.x() << ", " << projection.y() << ", " << projection.z() << " )\n";
0148 }
0149 }
0150 } else {
0151 if (!chamberGeometry)
0152 edm::LogWarning("MuonIdentification") << "Cannot get chamber geomtry for DetId: " << chamberId.rawId();
0153 if (!simUnitGeometry)
0154 edm::LogWarning("MuonIdentification") << "Cannot get detector unit geomtry for DetId: " << hit.detUnitId();
0155 }
0156 }
0157
0158 double MuonIdTruthInfo::matchChi2(const reco::Track& recoTrk, const SimTrack& simTrk) {
0159 double deltaPhi = fabs(recoTrk.phi() - simTrk.momentum().phi());
0160 if (deltaPhi > 1.8 * 3.1416)
0161 deltaPhi = 2 * 3.1416 - deltaPhi;
0162 return pow((recoTrk.p() - simTrk.momentum().rho()) / simTrk.momentum().rho(), 2) + pow(deltaPhi / (2 * 3.1416), 2) +
0163 pow((recoTrk.theta() - simTrk.momentum().theta()) / 3.1416, 2);
0164 }