CosmicMuonProducer/src/CosmicMuonLinksProducer.cc

0001 #include "RecoMuon/CosmicMuonProducer/src/CosmicMuonLinksProducer.h"
0002
0003 /**\class CosmicMuonLinksProducer
0004  *
0005  *
0006  * Original Author:  Chang Liu - Purdue University <chang.liu@cern.ch>
0007 **/
0008
0009 // system include files
0010 #include <memory>
0011
0012 // user include files
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024
0025 #include "DataFormats/TrackReco/interface/Track.h"
0026 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0027 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0028
0029 using namespace edm;
0030 using namespace std;
0031 using namespace reco;
0032
0033 //
0034 // constructors and destructor
0035 //
0036 CosmicMuonLinksProducer::CosmicMuonLinksProducer(const ParameterSet& iConfig) {
0037   category_ = "Muon|RecoMuon|CosmicMuon|CosmicMuonLinksProducer";
0038
0039   ParameterSet serviceParameters = iConfig.getParameter<ParameterSet>("ServiceParameters");
0040
0041   theService = new MuonServiceProxy(serviceParameters, consumesCollector());
0042
0043   std::vector<edm::ParameterSet> theMapPSets = iConfig.getParameter<std::vector<edm::ParameterSet> >("Maps");
0044   for (std::vector<edm::ParameterSet>::const_iterator iMPS = theMapPSets.begin(); iMPS != theMapPSets.end(); iMPS++) {
0045     edm::InputTag sTag = (*iMPS).getParameter<edm::InputTag>("subTrack");
0046     edm::InputTag pTag = (*iMPS).getParameter<edm::InputTag>("parentTrack");
0047
0048     edm::EDGetTokenT<reco::TrackCollection> subTrackTag = consumes<reco::TrackCollection>(sTag);
0049     edm::EDGetTokenT<reco::TrackCollection> parentTrackTag = consumes<reco::TrackCollection>(pTag);
0050     theTrackLinks.push_back(make_pair(subTrackTag, parentTrackTag));
0051     theTrackLinkNames.push_back(make_pair(sTag.label(), pTag.label()));
0052
0053     LogDebug(category_) << "preparing map between " << sTag << " & " << pTag;
0054     std::string mapname = sTag.label() + "To" + pTag.label();
0055     produces<reco::TrackToTrackMap>(mapname);
0056   }
0057 }
0058
0059 CosmicMuonLinksProducer::~CosmicMuonLinksProducer() {
0060   if (theService)
0061     delete theService;
0062 }
0063
0064 // ------------ method called to produce the data  ------------
0065 void CosmicMuonLinksProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0066   LogInfo(category_) << "Processing event number: " << iEvent.id();
0067
0068   theService->update(iSetup);
0069
0070   unsigned int counter =
0071       0;  ///DAMN I cannot read the label of the TOKEN so I need to do this stupid thing to create the labels of the products!
0072   for (std::vector<std::pair<edm::EDGetTokenT<reco::TrackCollection>,
0073                              edm::EDGetTokenT<reco::TrackCollection> > >::const_iterator iLink = theTrackLinks.begin();
0074        iLink != theTrackLinks.end();
0075        iLink++) {
0076 #ifdef EDM_ML_DEBUG
0077     edm::EDConsumerBase::Labels labels_first;
0078     edm::EDConsumerBase::Labels labels_second;
0079     labelsForToken(iLink->first, labels_first);
0080     labelsForToken(iLink->second, labels_second);
0081     LogDebug(category_) << "making map between " << labels_first.module << " and " << labels_second.module;
0082 #endif
0083     std::string mapname = theTrackLinkNames[counter].first + "To" + theTrackLinkNames[counter].second;
0084     reco::TrackToTrackMap ttmap;
0085
0086     Handle<reco::TrackCollection> subTracks;
0087     Handle<reco::TrackCollection> parentTracks;
0088
0089     iEvent.getByToken((*iLink).first, subTracks);
0090     iEvent.getByToken((*iLink).second, parentTracks);
0091
0092     ttmap = mapTracks(subTracks, parentTracks);
0093     LogTrace(category_) << "Mapped: " << theTrackLinkNames[counter].first << " " << subTracks->size() << " and "
0094                         << theTrackLinkNames[counter].second << " " << parentTracks->size()
0095                         << ", results: " << ttmap.size() << endl;
0096
0097     iEvent.put(std::make_unique<reco::TrackToTrackMap>(ttmap), mapname);
0098
0099     counter++;
0100   }
0101 }
0102
0103 reco::TrackToTrackMap CosmicMuonLinksProducer::mapTracks(const Handle<reco::TrackCollection>& subTracks,
0104                                                          const Handle<reco::TrackCollection>& parentTracks) const {
0105   reco::TrackToTrackMap map(subTracks, parentTracks);
0106   for (unsigned int position1 = 0; position1 != subTracks->size(); ++position1) {
0107     TrackRef track1(subTracks, position1);
0108     for (unsigned int position2 = 0; position2 != parentTracks->size(); ++position2) {
0109       TrackRef track2(parentTracks, position2);
0110       int shared = sharedHits(*track1, *track2);
0111       LogTrace(category_) << "sharedHits " << shared << " track1 " << track1->found() << " track2 " << track2->found()
0112                           << endl;
0113
0114       if (shared > (track1->found()) / 2)
0115         map.insert(track1, track2);
0116     }
0117   }
0118
0119   return map;
0120 }
0121
0122 int CosmicMuonLinksProducer::sharedHits(const reco::Track& track1, const reco::Track& track2) const {
0123   int match = 0;
0124
0125   for (trackingRecHit_iterator hit1 = track1.recHitsBegin(); hit1 != track1.recHitsEnd(); ++hit1) {
0126     if (!(*hit1)->isValid())
0127       continue;
0128     DetId id1 = (*hit1)->geographicalId();
0129     if (id1.det() != DetId::Muon)
0130       continue;  //ONLY MUON
0131     LogTrace(category_) << "first ID " << id1.rawId() << " " << (*hit1)->localPosition() << endl;
0132     GlobalPoint pos1 = theService->trackingGeometry()->idToDet(id1)->surface().toGlobal((*hit1)->localPosition());
0133
0134     for (trackingRecHit_iterator hit2 = track2.recHitsBegin(); hit2 != track2.recHitsEnd(); ++hit2) {
0135       if (!(*hit2)->isValid())
0136         continue;
0137
0138       DetId id2 = (*hit2)->geographicalId();
0139       if (id2.det() != DetId::Muon)
0140         continue;  //ONLY MUON
0141
0142       //          LogTrace(category_)<<"second ID "<<id2.rawId()<< (*hit2)->localPosition()<<endl;
0143
0144       if (id2.rawId() != id1.rawId())
0145         continue;
0146
0147       GlobalPoint pos2 = theService->trackingGeometry()->idToDet(id2)->surface().toGlobal((*hit2)->localPosition());
0148       if ((pos1 - pos2).mag() < 10e-5)
0149         match++;
0150     }
0151   }
0152
0153   return match;
0154 }