File indexing completed on 2023-03-17 11:17:16
0001 #include <memory>
0002 #include <iostream>
0003 #include <string>
0004 #include <map>
0005 #include <set>
0006
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "DataFormats/TrackReco/interface/Track.h"
0014 #include "DataFormats/MuonReco/interface/Muon.h"
0015 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0016 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0017 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0018
0019 namespace std {
0020 template <>
0021 class less<edm::RefToBase<reco::Track> > {
0022 public:
0023 bool operator()(const edm::RefToBase<reco::Track>& x, const edm::RefToBase<reco::Track>& y) const {
0024 return (x.id().id() < y.id().id()) || (x.key() < y.key()) || false;
0025 }
0026 };
0027 }
0028
0029 class testMuonAssociator : public edm::one::EDAnalyzer<> {
0030 public:
0031 explicit testMuonAssociator(const edm::ParameterSet& iConfig);
0032 virtual ~testMuonAssociator();
0033 virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
0034
0035 private:
0036 edm::EDGetTokenT<edm::View<reco::Track> > token_recoTracks;
0037 edm::EDGetTokenT<edm::View<reco::Track> > token_standAloneMuons;
0038 edm::EDGetTokenT<edm::View<reco::Track> > token_globalMuons;
0039 edm::EDGetTokenT<reco::MuonCollection> token_muons;
0040 edm::EDGetTokenT<TrackingParticleCollection> token_trackingTruth;
0041 edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> token_associatorByHits;
0042 edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> token_associatorByChi2;
0043 edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> token_associatorByPos;
0044 unsigned int m_flavour;
0045 double m_ptcut;
0046 };
0047
0048 std::ostream& operator<<(std::ostream& out, edm::RefToBase<reco::Track> ref) {
0049 out << std::fixed << " {" << std::setw(2) << ref->found() << "} "
0050 << " [" << std::setw(4) << ref.key() << "]"
0051 << " "
0052 << " pT: " << std::setw(6) << std::setprecision(3) << ref->pt() << " eta: " << std::setw(6)
0053 << std::setprecision(3) << ref->eta() << " phi: " << std::setw(6) << std::setprecision(3) << ref->phi();
0054 return out;
0055 }
0056
0057 std::ostream& operator<<(std::ostream& out, reco::MuonRef ref) {
0058 out << std::fixed;
0059 if (ref->isGlobalMuon()) {
0060 out << " {" << std::setw(2) << ref->innerTrack()->found() << "+" << std::setw(2) << ref->outerTrack()->found()
0061 << "} "
0062 << " [" << std::setw(4) << ref.key() << "]"
0063 << " "
0064 << " pT: " << std::setw(6) << std::setprecision(3) << ref->globalTrack()->pt() << " eta: " << std::setw(6)
0065 << std::setprecision(3) << ref->globalTrack()->eta() << " phi: " << std::setw(6) << std::setprecision(3)
0066 << ref->globalTrack()->phi();
0067 } else if (ref->isTrackerMuon()) {
0068 out << " {" << std::setw(2) << ref->innerTrack()->found() << " } "
0069 << " [" << std::setw(4) << ref.key() << "]"
0070 << " "
0071 << " pT: " << std::setw(6) << std::setprecision(3) << ref->innerTrack()->pt() << " eta: " << std::setw(6)
0072 << std::setprecision(3) << ref->innerTrack()->eta() << " phi: " << std::setw(6) << std::setprecision(3)
0073 << ref->innerTrack()->phi();
0074 } else if (ref->isStandAloneMuon()) {
0075 out << " { " << std::setw(2) << ref->outerTrack()->found() << "} "
0076 << " [" << std::setw(4) << ref.key() << "]"
0077 << " "
0078 << " pT: " << std::setw(6) << std::setprecision(3) << ref->outerTrack()->pt() << " eta: " << std::setw(6)
0079 << std::setprecision(3) << ref->outerTrack()->eta() << " phi: " << std::setw(6) << std::setprecision(3)
0080 << ref->outerTrack()->phi();
0081 } else {
0082 out << "(muon track not available)";
0083 }
0084 return out;
0085 }
0086
0087 std::ostream& operator<<(std::ostream& out, TrackingParticleRef ref) {
0088 out << std::fixed;
0089 out << " [" << std::setw(4) << ref.key() << "]"
0090 << " type:" << std::setw(6) << ref->pdgId() << " pT: " << std::setw(6) << std::setprecision(3) << ref->pt()
0091 << " eta: " << std::setw(6) << std::setprecision(3) << ref->eta() << " phi: " << std::setw(6)
0092 << std::setprecision(3) << ref->phi();
0093 return out;
0094 }
0095
0096 struct Quality {
0097 double byChi2;
0098 double byHits;
0099 double byPosition;
0100 };
0101
0102 template <class Ref>
0103 class Associations {
0104 public:
0105 typedef std::map<Ref, Quality> map_type;
0106 typedef std::vector<std::pair<Ref, double> > association_type;
0107
0108 template <class Key, class Associator>
0109 Associations(const std::string& label,
0110 const Key& candidate,
0111 const Associator& byHits,
0112 const Associator& byChi2,
0113 const Associator& byPosition)
0114 : m_map(), m_label(label) {
0115 if (byHits.find(candidate) != byHits.end())
0116 fillByHits(byHits[candidate]);
0117 if (byChi2.find(candidate) != byChi2.end())
0118 fillByChi2(byChi2[candidate]);
0119 if (byPosition.find(candidate) != byPosition.end())
0120 fillByPosition(byPosition[candidate]);
0121 }
0122
0123 void fillByHits(const association_type& found) {
0124 for (typename association_type::const_iterator it = found.begin(); it != found.end(); ++it) {
0125 const Ref& ref = it->first;
0126 m_map.insert(std::make_pair(ref, Quality()));
0127 m_map[ref].byHits = it->second;
0128 }
0129 }
0130
0131 void fillByChi2(const association_type& found) {
0132 for (typename association_type::const_iterator it = found.begin(); it != found.end(); ++it) {
0133 const Ref& ref = it->first;
0134 m_map.insert(std::make_pair(ref, Quality()));
0135 m_map[ref].byChi2 = -it->second;
0136 }
0137 }
0138
0139 void fillByPosition(const association_type& found) {
0140 for (typename association_type::const_iterator it = found.begin(); it != found.end(); ++it) {
0141 const Ref& ref = it->first;
0142 m_map.insert(std::make_pair(ref, Quality()));
0143 m_map[ref].byPosition = -it->second;
0144 }
0145 }
0146
0147 const map_type& map(void) const { return m_map; }
0148
0149 void dump(std::ostream& out) const {
0150 std::stringstream buffer;
0151 for (typename map_type::const_iterator it = m_map.begin(); it != m_map.end(); ++it) {
0152 buffer << " " << std::setw(7) << std::left << m_label << std::right << it->first;
0153 if (it->second.byHits)
0154 buffer << " [" << std::setw(6) << std::setprecision(3) << it->second.byHits << "]";
0155 else
0156 buffer << " ";
0157 if (it->second.byChi2)
0158 buffer << " [" << std::setw(6) << std::setprecision(3) << it->second.byChi2 << "]";
0159 else
0160 buffer << " ";
0161 if (it->second.byPosition)
0162 buffer << " [" << std::setw(6) << std::setprecision(3) << it->second.byPosition << "]";
0163 else
0164 buffer << " ";
0165 buffer << std::endl;
0166 }
0167 out << buffer.str();
0168 }
0169
0170 private:
0171 map_type m_map;
0172 std::string m_label;
0173 };
0174
0175 template <class Ref>
0176 std::ostream& operator<<(std::ostream& out, const Associations<Ref>& association) {
0177 association.dump(out);
0178 return out;
0179 }
0180
0181 template <typename T>
0182 std::ostream& operator<<(std::ostream& out, const std::pair<edm::Ref<T>, double>& assoc) {
0183 out << assoc.first << " quality: " << assoc.second;
0184 return out;
0185 }
0186
0187 testMuonAssociator::testMuonAssociator(edm::ParameterSet const& iConfig) {
0188 token_recoTracks = consumes<edm::View<reco::Track> >(iConfig.getParameter<edm::InputTag>("tracks"));
0189 token_standAloneMuons =
0190 consumes<edm::View<reco::Track> >(iConfig.getParameter<edm::InputTag>("standAloneMuonTracks"));
0191 token_globalMuons = consumes<edm::View<reco::Track> >(iConfig.getParameter<edm::InputTag>("globalMuonTracks"));
0192 token_muons = consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons"));
0193 token_trackingTruth = consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingTruth"));
0194 token_associatorByHits = consumes<reco::TrackToTrackingParticleAssociator>(edm::InputTag("trackAssociatorByHits"));
0195 token_associatorByChi2 = consumes<reco::TrackToTrackingParticleAssociator>(edm::InputTag("trackAssociatorByChi2"));
0196 token_associatorByPos = consumes<reco::TrackToTrackingParticleAssociator>(edm::InputTag("trackAssociatorByPosition"));
0197 m_flavour = iConfig.getParameter<unsigned int>("leptonFlavour");
0198 m_ptcut = iConfig.getParameter<double>("minPt");
0199 }
0200
0201 testMuonAssociator::~testMuonAssociator() {}
0202
0203 void testMuonAssociator::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0204
0205 edm::Handle<reco::TrackToTrackingParticleAssociator> associatorByHits;
0206 event.getByToken(token_associatorByHits, associatorByHits);
0207
0208 edm::Handle<reco::TrackToTrackingParticleAssociator> associatorByChi2;
0209 event.getByToken(token_associatorByChi2, associatorByChi2);
0210
0211 edm::Handle<reco::TrackToTrackingParticleAssociator> associatorByPos;
0212 event.getByToken(token_associatorByPos, associatorByPos);
0213
0214
0215 edm::Handle<edm::View<reco::Track> > recoTrackHandle;
0216 event.getByToken(token_recoTracks, recoTrackHandle);
0217 const edm::View<reco::Track>& recoTrackCollection = *(recoTrackHandle.product());
0218
0219 edm::Handle<edm::View<reco::Track> > standAloneMuonHandle;
0220 event.getByToken(token_standAloneMuons, standAloneMuonHandle);
0221 const edm::View<reco::Track>& standAloneMuonCollection = *(standAloneMuonHandle.product());
0222
0223 edm::Handle<edm::View<reco::Track> > globalMuonTrackHandle;
0224 event.getByToken(token_globalMuons, globalMuonTrackHandle);
0225 const edm::View<reco::Track>& globalMuonTrackCollection = *(globalMuonTrackHandle.product());
0226
0227 edm::Handle<reco::MuonCollection> muonHandle;
0228 event.getByToken(token_muons, muonHandle);
0229 const reco::MuonCollection& muonCollection = *(muonHandle.product());
0230
0231 edm::Handle<TrackingParticleCollection> trackingParticleHandle;
0232 event.getByToken(token_trackingTruth, trackingParticleHandle);
0233 const TrackingParticleCollection& trackingParticleCollection = *(trackingParticleHandle.product());
0234
0235 std::cout << std::fixed;
0236
0237 std::cout << std::endl;
0238 std::cout << "Found " << std::setw(6) << trackingParticleCollection.size() << " TrackingParticles" << std::flush;
0239 unsigned int count = 0;
0240 for (TrackingParticleCollection::size_type i = 0; i < trackingParticleCollection.size(); ++i)
0241 if ((std::abs(trackingParticleCollection[i].pdgId()) == (int)m_flavour) and
0242 (trackingParticleCollection[i].pt() >= m_ptcut))
0243 ++count;
0244
0245 std::cout << " ( " << std::setw(2) << count << " leptons with pT above " << m_ptcut << " GeV)" << std::endl;
0246 std::cout << " " << std::setw(6) << recoTrackCollection.size() << " Tracks" << std::endl;
0247 std::cout << " " << std::setw(6) << muonCollection.size() << " Global muons" << std::endl;
0248 std::cout << " " << std::setw(6) << standAloneMuonCollection.size() << " StandAlone muons" << std::endl;
0249 std::cout << std::endl;
0250
0251
0252 {
0253 reco::SimToRecoCollection bychi2_tracks;
0254 reco::SimToRecoCollection bychi2_globaltrack;
0255 reco::SimToRecoCollection bychi2_standalone;
0256 reco::SimToRecoCollection bypos_tracks;
0257 reco::SimToRecoCollection bypos_globaltrack;
0258 reco::SimToRecoCollection bypos_standalone;
0259 reco::SimToRecoCollection byhits_tracks;
0260 reco::SimToRecoCollection byhits_globaltrack;
0261 reco::SimToRecoCollection byhits_standalone;
0262 try {
0263 bychi2_tracks = associatorByChi2->associateSimToReco(recoTrackHandle, trackingParticleHandle);
0264 } catch (const std::exception& e) {
0265 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0266 }
0267 try {
0268 bychi2_globaltrack = associatorByChi2->associateSimToReco(globalMuonTrackHandle, trackingParticleHandle);
0269 } catch (const std::exception& e) {
0270 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0271 }
0272 try {
0273 bychi2_standalone = associatorByChi2->associateSimToReco(standAloneMuonHandle, trackingParticleHandle);
0274 } catch (const std::exception& e) {
0275 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0276 }
0277 try {
0278 byhits_tracks = associatorByHits->associateSimToReco(recoTrackHandle, trackingParticleHandle);
0279 } catch (const std::exception& e) {
0280 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0281 }
0282 try {
0283 byhits_globaltrack = associatorByHits->associateSimToReco(globalMuonTrackHandle, trackingParticleHandle);
0284 } catch (const std::exception& e) {
0285 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0286 }
0287 try {
0288 byhits_standalone = associatorByHits->associateSimToReco(standAloneMuonHandle, trackingParticleHandle);
0289 } catch (const std::exception& e) {
0290 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0291 }
0292
0293
0294
0295
0296
0297
0298
0299
0300
0301
0302
0303
0304 for (TrackingParticleCollection::size_type i = 0; i < trackingParticleCollection.size(); ++i) {
0305 TrackingParticleRef tp(trackingParticleHandle, i);
0306 if ((std::abs(tp->pdgId()) != (int)m_flavour) or (tp->pt() < m_ptcut))
0307 continue;
0308 std::cout << "--> TrackingParticle" << tp << std::endl;
0309 std::cout << Associations<edm::RefToBase<reco::Track> >("Track", tp, byhits_tracks, bychi2_tracks, bypos_tracks);
0310 std::cout << Associations<edm::RefToBase<reco::Track> >(
0311 "Local", tp, byhits_standalone, bychi2_standalone, bypos_standalone);
0312 std::cout << Associations<edm::RefToBase<reco::Track> >(
0313 "Global", tp, byhits_globaltrack, bychi2_globaltrack, bypos_globaltrack);
0314 }
0315 }
0316
0317
0318 reco::RecoToSimCollection byhits_muon;
0319 reco::RecoToSimCollection bychi2_muon;
0320 reco::RecoToSimCollection bypos_muon;
0321 try {
0322 byhits_muon = associatorByHits->associateRecoToSim(recoTrackHandle, trackingParticleHandle);
0323 } catch (const std::exception& e) {
0324 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0325 }
0326 try {
0327 bychi2_muon = associatorByChi2->associateRecoToSim(recoTrackHandle, trackingParticleHandle);
0328 } catch (const std::exception& e) {
0329 std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl;
0330 }
0331
0332
0333
0334
0335
0336 for (reco::MuonCollection::size_type i = 0; i < muonCollection.size(); ++i) {
0337 reco::MuonRef muon(muonHandle, i);
0338 std::cout << "<-- Muon " << muon << std::endl;
0339 std::cout << Associations<TrackingParticleRef>(
0340 "TrackingParticle", edm::RefToBase<reco::Track>(muon->track()), byhits_muon, bychi2_muon, bypos_muon);
0341 }
0342
0343
0344 reco::RecoToSimCollection byhits_global;
0345 reco::RecoToSimCollection bychi2_global;
0346 reco::RecoToSimCollection bypos_global;
0347 try {
0348 byhits_global = associatorByHits->associateRecoToSim(globalMuonTrackHandle, trackingParticleHandle);
0349 } catch (const std::exception& e) {
0350 std::cerr << std::endl << "Error building reco::RecoToSimCollection:\n" << e.what() << std::endl;
0351 }
0352 try {
0353 bychi2_global = associatorByChi2->associateRecoToSim(globalMuonTrackHandle, trackingParticleHandle);
0354 } catch (const std::exception& e) {
0355 std::cerr << std::endl << "Error building reco::RecoToSimCollection:\n" << e.what() << std::endl;
0356 }
0357
0358
0359
0360
0361
0362 for (edm::View<reco::Track>::size_type i = 0; i < globalMuonTrackCollection.size(); ++i) {
0363 edm::RefToBase<reco::Track> track(globalMuonTrackHandle, i);
0364 std::cout << "<-- Local " << track << std::endl;
0365 std::cout << Associations<TrackingParticleRef>(
0366 "TrackingParticle", track, byhits_global, bychi2_global, bypos_global);
0367 }
0368
0369
0370 reco::RecoToSimCollection byhits_standalone;
0371 reco::RecoToSimCollection bychi2_standalone;
0372 reco::RecoToSimCollection bypos_standalone;
0373 try {
0374 byhits_standalone = associatorByHits->associateRecoToSim(standAloneMuonHandle, trackingParticleHandle);
0375 } catch (const std::exception& e) {
0376 std::cerr << std::endl << "Error building reco::RecoToSimCollection:\n" << e.what() << std::endl;
0377 }
0378 try {
0379 bychi2_standalone = associatorByChi2->associateRecoToSim(standAloneMuonHandle, trackingParticleHandle);
0380 } catch (const std::exception& e) {
0381 std::cerr << std::endl << "Error building reco::RecoToSimCollection:\n" << e.what() << std::endl;
0382 }
0383
0384
0385
0386
0387
0388 for (edm::View<reco::Track>::size_type i = 0; i < standAloneMuonCollection.size(); ++i) {
0389 edm::RefToBase<reco::Track> track(standAloneMuonHandle, i);
0390 std::cout << "<-- Local " << track << std::endl;
0391 std::cout << Associations<TrackingParticleRef>(
0392 "TrackingParticle", track, byhits_standalone, bychi2_standalone, bypos_standalone);
0393 }
0394
0395 std::cout << std::endl;
0396 }
0397
0398 #include "FWCore/Framework/interface/MakerMacros.h"
0399 DEFINE_FWK_MODULE(testMuonAssociator);