/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoMuon/MuonIdentification/src/MuonIdTruthInfo.cc

Go to the documentation of this file.

#include "RecoMuon/MuonIdentification/interface/MuonIdTruthInfo.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/DTChamberId.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"

void MuonIdTruthInfo::truthMatchMuon(const edm::Event& iEvent, 
                                     const edm::EventSetup& iSetup,
                                     reco::Muon& aMuon)
{
   // get a list of simulated track and find a track with the best match to
   // the muon.track(). Use its id and chamber id to localize hits
   // If a hit has non-zero local z coordinate, it's position wrt
   // to the center of a chamber is extrapolated by a straight line
   
   edm::Handle<edm::SimTrackContainer> simTracks;
   iEvent.getByType<edm::SimTrackContainer>(simTracks);
   if (! simTracks.isValid() ) {
      LogTrace("MuonIdentification") <<"No tracks found";
      return;
   }
   
   // get the tracking Geometry
   edm::ESHandle<GlobalTrackingGeometry> geometry;
   iSetup.get<GlobalTrackingGeometryRecord>().get(geometry);
   if ( ! geometry.isValid() )
     throw cms::Exception("FatalError") << "Unable to find GlobalTrackingGeometryRecord in event!\n";
   
   float bestMatchChi2 = 9999; //minimization creteria
   unsigned int bestMatch = 0;
   const unsigned int offset = 0; // kludge to fix a problem in trackId matching between tracks and hits.
   
   for( edm::SimTrackContainer::const_iterator simTrk = simTracks->begin();
        simTrk != simTracks->end(); simTrk++ )
     {
        float chi2 = matchChi2(*aMuon.track().get(),*simTrk);
        if (chi2>bestMatchChi2) continue;
        bestMatchChi2 = chi2;
        bestMatch = simTrk->trackId();
     }
   
   bestMatch -= offset;
   
   std::vector<reco::MuonChamberMatch>& matches = aMuon.matches();
   int numberOfTruthMatchedChambers = 0;

   // loop over chambers
   for(std::vector<reco::MuonChamberMatch>::iterator chamberMatch = matches.begin();
       chamberMatch != matches.end(); chamberMatch ++)
     {
        if (chamberMatch->id.det() != DetId::Muon ) {
           edm::LogWarning("MuonIdentification") << "Detector id of a muon chamber corresponds to not a muon detector";
           continue;
        }
        reco::MuonSegmentMatch bestSegmentMatch;
        double distance = 99999;
        
        if ( chamberMatch->id.subdetId() == MuonSubdetId::DT) {
           DTChamberId detId(chamberMatch->id.rawId());

           edm::Handle<edm::PSimHitContainer> simHits;
           iEvent.getByLabel("g4SimHits", "MuonDTHits", simHits);
           if ( simHits.isValid() ) {
              for( edm::PSimHitContainer::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++)
                if ( hit->trackId() == bestMatch ) checkSimHitForBestMatch(bestSegmentMatch, distance, *hit, detId, geometry );
           }else LogTrace("MuonIdentification") <<"No DT simulated hits are found";
        }

        if ( chamberMatch->id.subdetId() == MuonSubdetId::CSC) {
           CSCDetId detId(chamberMatch->id.rawId());

           edm::Handle<edm::PSimHitContainer> simHits;
           iEvent.getByLabel("g4SimHits", "MuonCSCHits", simHits);
           if ( simHits.isValid() ) {
              for( edm::PSimHitContainer::const_iterator hit = simHits->begin(); hit != simHits->end(); hit++)
                if ( hit->trackId() == bestMatch ) checkSimHitForBestMatch(bestSegmentMatch, distance, *hit, detId, geometry );
           }else LogTrace("MuonIdentification") <<"No CSC simulated hits are found";
        }
        if (distance < 9999) {
           chamberMatch->truthMatches.push_back( bestSegmentMatch );
           numberOfTruthMatchedChambers++;
           LogTrace("MuonIdentification") << "Best truth matched hit:" <<
              "\tDetId: " << chamberMatch->id.rawId() << "\n" <<
              "\tprojection: ( " << bestSegmentMatch.x << ", " << bestSegmentMatch.y << " )\n";
        }
     }
   LogTrace("MuonIdentification") << "Truth matching summary:\n\tnumber of chambers: " << matches.size() <<
     "\n\tnumber of truth matched chambers: " << numberOfTruthMatchedChambers << "\n";
}

void MuonIdTruthInfo::checkSimHitForBestMatch(reco::MuonSegmentMatch& segmentMatch,
                                              double& distance,
                                              const PSimHit& hit, 
                                              const DetId& chamberId,
                                              const edm::ESHandle<GlobalTrackingGeometry>& geometry)
{
   // find the hit position projection at the reference surface of the chamber:
   // first get entry and exit point of the hit in the global coordinates, then
   // get local coordinates of these points wrt the chamber and then find the
   // projected X-Y coordinates
   
   const GeomDet* chamberGeometry = geometry->idToDet( chamberId );
   const GeomDet* simUnitGeometry = geometry->idToDet( DetId(hit.detUnitId()) );
   
   if (chamberGeometry && simUnitGeometry ) {
      LocalPoint entryPoint = chamberGeometry->toLocal( simUnitGeometry->toGlobal( hit.entryPoint() ) );
      LocalPoint exitPoint =  chamberGeometry->toLocal( simUnitGeometry->toGlobal( hit.exitPoint()  ) );
      LocalVector direction = exitPoint - entryPoint;
      if ( fabs(direction.z()) > 0.001) {
         LocalPoint projection = entryPoint - direction*(entryPoint.z()/direction.z());
         if ( fabs(projection.z()) > 0.001 ) 
           edm::LogWarning("MuonIdentification") << "z coordinate of the hit projection must be zero and it's not!\n";
         
         double new_distance = 99999;
         if( entryPoint.z()*exitPoint.z() < -1 ) // changed sign, so the reference point is inside
           new_distance = 0;
         else {
            if ( fabs(entryPoint.z()) < fabs(exitPoint.z()) )
              new_distance = fabs(entryPoint.z());
            else
              new_distance = fabs(exitPoint.z());
         }
         
         if (new_distance < distance) { 
            // find a SimHit closer to the reference surface, update segmentMatch
            segmentMatch.x = projection.x();
            segmentMatch.y = projection.y();
            segmentMatch.xErr = 0;
            segmentMatch.yErr = 0;
            segmentMatch.dXdZ = direction.x()/direction.z();
            segmentMatch.dYdZ = direction.y()/direction.z();
            segmentMatch.dXdZErr = 0;
            segmentMatch.dYdZErr = 0;
            distance = new_distance;
            LogTrace("MuonIdentificationVerbose") << "Better truth matched segment found:\n" <<
              "\tDetId: " << chamberId.rawId() << "\n" <<
              "\tentry point: ( " << entryPoint.x() << ", " << entryPoint.y() << ", " << entryPoint.z() << " )\n" <<
              "\texit point: ( " << exitPoint.x() << ", " << exitPoint.y() << ", " << exitPoint.z() << " )\n" <<
              "\tprojection: ( " << projection.x() << ", " << projection.y() << ", " << projection.z() << " )\n";
         }
      }
   }else{
      if ( ! chamberGeometry ) edm::LogWarning("MuonIdentification") << "Cannot get chamber geomtry for DetId: " << chamberId.rawId();
      if ( ! simUnitGeometry ) edm::LogWarning("MuonIdentification") << "Cannot get detector unit geomtry for DetId: " << hit.detUnitId();
   }
}

double MuonIdTruthInfo::matchChi2( const reco::Track& recoTrk, const SimTrack& simTrk)
{
   double deltaPhi = fabs(recoTrk.phi()-simTrk.momentum().phi());
   if (deltaPhi>1.8*3.1416) deltaPhi = 2*3.1416-deltaPhi; // take care of phi discontinuity
   return pow((recoTrk.p()-simTrk.momentum().rho())/simTrk.momentum().rho(),2)+
     pow(deltaPhi/(2*3.1416),2)+pow((recoTrk.theta()-simTrk.momentum().theta())/3.1416,2);
}