TrackerMuonHitExtractor.cc

Go to the documentation of this file.

//
// modified & integrated by Giovanni Abbiendi
// from code by Arun Luthra: UserCode/luthra/MuonTrackSelector/src/MuonTrackSelector.cc
//
#include "SimMuon/MCTruth/interface/TrackerMuonHitExtractor.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/DTChamberId.h"
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
#include "DataFormats/MuonReco/interface/MuonSelectors.h"
#include <sstream>

TrackerMuonHitExtractor::TrackerMuonHitExtractor(const edm::ParameterSet& parset,
                         edm::ConsumesCollector && ic) :
  inputDTRecSegment4DToken_(ic.consumes<DTRecSegment4DCollection>(parset.getParameter<edm::InputTag>("inputDTRecSegment4DCollection"))),
  inputCSCSegmentToken_(ic.consumes<CSCSegmentCollection>(parset.getParameter<edm::InputTag>("inputCSCSegmentCollection"))),
  inputDTRecSegment4DCollection_(parset.getParameter<edm::InputTag>("inputDTRecSegment4DCollection")),
  inputCSCSegmentCollection_(parset.getParameter<edm::InputTag>("inputCSCSegmentCollection"))
{
}

TrackerMuonHitExtractor::TrackerMuonHitExtractor(const edm::ParameterSet& parset ) :
  inputDTRecSegment4DCollection_(parset.getParameter<edm::InputTag>("inputDTRecSegment4DCollection")),
  inputCSCSegmentCollection_(parset.getParameter<edm::InputTag>("inputCSCSegmentCollection"))
{
}

TrackerMuonHitExtractor::~TrackerMuonHitExtractor() {
}

void TrackerMuonHitExtractor::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) 
{
  iEvent.getByLabel(inputDTRecSegment4DCollection_, dtSegmentCollectionH_);
  iEvent.getByLabel(inputCSCSegmentCollection_, cscSegmentCollectionH_);

  edm::LogVerbatim("TrackerMuonHitExtractor") <<"\nThere are "<< dtSegmentCollectionH_->size()<<" DT segments.";
  unsigned int index_dt_segment = 0;
  for(DTRecSegment4DCollection::const_iterator segment = dtSegmentCollectionH_->begin();
      segment != dtSegmentCollectionH_->end(); ++segment , index_dt_segment++) {
    LocalPoint  segmentLocalPosition       = segment->localPosition();
    LocalVector segmentLocalDirection      = segment->localDirection();
    LocalError  segmentLocalPositionError  = segment->localPositionError();
    LocalError  segmentLocalDirectionError = segment->localDirectionError();
    DetId geoid = segment->geographicalId();
    DTChamberId dtdetid = DTChamberId(geoid);
    int wheel = dtdetid.wheel();
    int station = dtdetid.station();
    int sector = dtdetid.sector();
    
    float segmentX = segmentLocalPosition.x();
    float segmentY = segmentLocalPosition.y();
    float segmentdXdZ = segmentLocalDirection.x()/segmentLocalDirection.z();
    float segmentdYdZ = segmentLocalDirection.y()/segmentLocalDirection.z();
    float segmentXerr = sqrt(segmentLocalPositionError.xx());
    float segmentYerr = sqrt(segmentLocalPositionError.yy());
    float segmentdXdZerr = sqrt(segmentLocalDirectionError.xx());
    float segmentdYdZerr = sqrt(segmentLocalDirectionError.yy());

    edm::LogVerbatim("TrackerMuonHitExtractor") 
      <<"\nDT segment index :"<<index_dt_segment
      <<"\nchamber Wh:"<<wheel<<",St:"<<station<<",Se:"<<sector
      <<"\nLocal Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
      <<"Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
  }

  edm::LogVerbatim("TrackerMuonHitExtractor") <<"\nThere are "<< cscSegmentCollectionH_->size()<<" CSC segments.";
  unsigned int index_csc_segment = 0;
  for(CSCSegmentCollection::const_iterator segment = cscSegmentCollectionH_->begin();
      segment != cscSegmentCollectionH_->end(); ++segment , index_csc_segment++) {
    LocalPoint  segmentLocalPosition       = segment->localPosition();
    LocalVector segmentLocalDirection      = segment->localDirection();
    LocalError  segmentLocalPositionError  = segment->localPositionError();
    LocalError  segmentLocalDirectionError = segment->localDirectionError();

    DetId geoid = segment->geographicalId();
    CSCDetId cscdetid = CSCDetId(geoid);
    int endcap = cscdetid.endcap();
    int station = cscdetid.station();
    int ring = cscdetid.ring();
    int chamber = cscdetid.chamber(); 
    
    float segmentX = segmentLocalPosition.x();
    float segmentY = segmentLocalPosition.y();
    float segmentdXdZ = segmentLocalDirection.x()/segmentLocalDirection.z();
    float segmentdYdZ = segmentLocalDirection.y()/segmentLocalDirection.z();
    float segmentXerr = sqrt(segmentLocalPositionError.xx());
    float segmentYerr = sqrt(segmentLocalPositionError.yy());
    float segmentdXdZerr = sqrt(segmentLocalDirectionError.xx());
    float segmentdYdZerr = sqrt(segmentLocalDirectionError.yy());

    edm::LogVerbatim("TrackerMuonHitExtractor") 
      <<"\nCSC segment index :"<<index_csc_segment
      <<"\nchamber Endcap:"<<endcap<<",St:"<<station<<",Ri:"<<ring<<",Ch:"<<chamber
      <<"\nLocal Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
      <<"Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
  }

}
std::vector<const TrackingRecHit *>
TrackerMuonHitExtractor::getMuonHits(const reco::Muon &mu) const {
        std::vector<const TrackingRecHit *> ret;

    int wheel, station, sector;
    int endcap, /*station, */ ring, chamber;
    
    edm::LogVerbatim("TrackerMuonHitExtractor") <<"Number of chambers: "<<mu.matches().size()
                          <<", arbitrated: "<<mu.numberOfMatches(reco::Muon::SegmentAndTrackArbitration)
                          <<", tot (no arbitration): "<<mu.numberOfMatches(reco::Muon::NoArbitration);
    unsigned int index_chamber = 0;
    int n_segments_noArb = 0;
    
    for(std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = mu.matches().begin();
        chamberMatch != mu.matches().end(); ++chamberMatch, index_chamber++) {
      std::stringstream chamberStr;
      chamberStr <<"\nchamber index: "<<index_chamber; 
      
      int subdet = chamberMatch->detector();
      DetId did = chamberMatch->id;
      
      if (subdet == MuonSubdetId::DT) {
        DTChamberId dtdetid = DTChamberId(did);
        wheel = dtdetid.wheel();
        station = dtdetid.station();
        sector = dtdetid.sector();
        chamberStr << ", DT chamber Wh:"<<wheel<<",St:"<<station<<",Se:"<<sector;
      } 
      else if (subdet == MuonSubdetId::CSC) {
        CSCDetId cscdetid = CSCDetId(did);
        endcap = cscdetid.endcap();
        station = cscdetid.station();
        ring = cscdetid.ring();
        chamber = cscdetid.chamber();
        chamberStr << ", CSC chamber End:"<<endcap<<",St:"<<station<<",Ri:"<<ring<<",Ch:"<<chamber;
      }
      
      chamberStr << ", Number of segments: "<<chamberMatch->segmentMatches.size();
      edm::LogVerbatim("TrackerMuonHitExtractor") << chamberStr.str();
      n_segments_noArb = n_segments_noArb + chamberMatch->segmentMatches.size();

      unsigned int index_segment = 0;
      
      for(std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();
          segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch, index_segment++) {
        
        float segmentX = segmentMatch->x;
        float segmentY = segmentMatch->y ;
        float segmentdXdZ = segmentMatch->dXdZ;
        float segmentdYdZ = segmentMatch->dYdZ;
        float segmentXerr = segmentMatch->xErr;
        float segmentYerr = segmentMatch->yErr;
        float segmentdXdZerr = segmentMatch->dXdZErr;
        float segmentdYdZerr = segmentMatch->dYdZErr;
        
        CSCSegmentRef segmentCSC = segmentMatch->cscSegmentRef;
        DTRecSegment4DRef segmentDT = segmentMatch->dtSegmentRef;
        
        bool segment_arbitrated_Ok = (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) && 
                      segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR));
        
        std::string ARBITRATED(" ***Arbitrated Off*** ");
        if (segment_arbitrated_Ok) ARBITRATED = " ***ARBITRATED OK*** ";

        if (subdet == MuonSubdetId::DT) {         
          edm::LogVerbatim("TrackerMuonHitExtractor")
        <<"\n\t segment index: "<<index_segment << ARBITRATED
        <<"\n\t  Local Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
        <<"\n\t  Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
          
          // with the following line the DT segments failing standard arbitration are skipped 
          if (!segment_arbitrated_Ok) continue;

              if (segmentDT.get() != nullptr) {
        const DTRecSegment4D* segment = segmentDT.get();
        
        edm::LogVerbatim("TrackerMuonHitExtractor")<<"\t ===> MATCHING with DT segment with index = "<<segmentDT.key();
        
        if(segment->hasPhi()) {
          const DTChamberRecSegment2D* phiSeg = segment->phiSegment();
          std::vector<const TrackingRecHit*> phiHits = phiSeg->recHits();
          for(std::vector<const TrackingRecHit*>::const_iterator ihit = phiHits.begin();
              ihit != phiHits.end(); ++ihit) {
                    ret.push_back(*ihit);
          }
        }
        
        if(segment->hasZed()) {
          const DTSLRecSegment2D* zSeg = (*segment).zSegment();
          std::vector<const TrackingRecHit*> zedHits = zSeg->recHits();
          for(std::vector<const TrackingRecHit*>::const_iterator ihit = zedHits.begin();
              ihit != zedHits.end(); ++ihit) {
                    ret.push_back(*ihit);
          }
        }
          } else edm::LogWarning("TrackerMuonHitExtractor")<<"\n***WARNING: UNMATCHED DT segment ! \n";
            } // if (subdet == MuonSubdetId::DT)
        
        else if (subdet == MuonSubdetId::CSC) {
          edm::LogVerbatim("TrackerMuonHitExtractor")
        <<"\n\t segment index: "<<index_segment << ARBITRATED
        <<"\n\t  Local Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
        <<"\n\t  Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
          
          // with the following line the CSC segments failing standard arbitration are skipped 
          if (!segment_arbitrated_Ok) continue;

              if (segmentCSC.get() != nullptr) {
        const CSCSegment* segment = segmentCSC.get();
        
        edm::LogVerbatim("TrackerMuonHitExtractor")<<"\t ===> MATCHING with CSC segment with index = "<<segmentCSC.key();
        
        std::vector<const TrackingRecHit*> hits = segment->recHits();
        for(std::vector<const TrackingRecHit*>::const_iterator ihit = hits.begin();
            ihit != hits.end(); ++ihit) {
                  ret.push_back(*ihit);
        }
          } else edm::LogWarning("TrackerMuonHitExtractor")<<"\n***WARNING: UNMATCHED CSC segment ! \n";
        }   // else if (subdet == MuonSubdetId::CSC)
        
      } // loop on vector<MuonSegmentMatch>   
    }  // loop on vector<MuonChamberMatch>  

    edm::LogVerbatim("TrackerMuonHitExtractor")<<"\n N. matched Segments before arbitration = "<<n_segments_noArb;

        return ret;
}