MuonIdentification/src/MuonIdTruthInfo.cc

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   // get a list of simulated track and find a track with the best match to
0021   // the muon.track(). Use its id and chamber id to localize hits
0022   // If a hit has non-zero local z coordinate, it's position wrt
0023   // to the center of a chamber is extrapolated by a straight line
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;  //minimization creteria
0033   unsigned int bestMatch = 0;
0034   const unsigned int offset = 0;  // kludge to fix a problem in trackId matching between tracks and hits.
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   // loop over chambers
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   // find the hit position projection at the reference surface of the chamber:
0105   // first get entry and exit point of the hit in the global coordinates, then
0106   // get local coordinates of these points wrt the chamber and then find the
0107   // projected X-Y coordinates
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)  // changed sign, so the reference point is inside
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         // find a SimHit closer to the reference surface, update segmentMatch
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;  // take care of phi discontinuity
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 }