File indexing completed on 2024-04-06 12:24:36
0001 #include <iostream>
0002 #include <map>
0003 #include <memory>
0004 #include <string>
0005 #include <set>
0006 #include <tuple>
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/MuonReco/interface/Muon.h"
0016 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0018 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0019
0020 namespace std {
0021 template <>
0022 class less<edm::RefToBase<reco::Track> > {
0023 public:
0024 bool operator()(const edm::RefToBase<reco::Track>& x, const edm::RefToBase<reco::Track>& y) const {
0025 return (x.id().processIndex() < y.id().processIndex()) || (x.id().productIndex() < y.id().productIndex()) ||
0026 (x.key() < y.key()) || false;
0027 }
0028 };
0029 }
0030
0031 class testLeptonAssociator : public edm::one::EDAnalyzer<> {
0032 public:
0033 explicit testLeptonAssociator(const edm::ParameterSet& iConfig);
0034 virtual void analyze(const edm::Event& iEvent, const edm::EventSetup& setup) override;
0035
0036 private:
0037 edm::EDGetTokenT<edm::View<reco::Track> > token_recoTracks;
0038 edm::EDGetTokenT<edm::View<reco::Track> > token_standAloneMuons;
0039 edm::EDGetTokenT<edm::View<reco::Track> > token_globalMuons;
0040 edm::EDGetTokenT<reco::MuonCollection> token_muons;
0041 edm::EDGetTokenT<TrackingParticleCollection> token_trackingTruth;
0042 edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> token_associatorByHits;
0043 edm::EDGetTokenT<reco::TrackToTrackingParticleAssociator> token_associatorByChi2;
0044
0045 unsigned int m_flavour;
0046 double m_ptcut;
0047 };
0048
0049 std::ostream& operator<<(std::ostream& out, edm::RefToBase<reco::Track> ref) {
0050 out << " {" << std::setw(2) << ref->found() << "} "
0051 << " [" << std::setw(4) << ref.key() << "]"
0052 << " "
0053 << " pT: " << std::setw(6) << std::setprecision(3) << ref->pt() << " eta: " << std::setw(6)
0054 << std::setprecision(3) << ref->eta() << " phi: " << std::setw(6) << std::setprecision(3) << ref->phi();
0055 return out;
0056 }
0057
0058 std::ostream& operator<<(std::ostream& out, reco::MuonRef ref) {
0059 if (ref->isGlobalMuon()) {
0060 out << " {" << std::setw(2) << ref->track()->found() << "+" << std::setw(2) << ref->standAloneMuon()->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->track()->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->standAloneMuon()->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 << "(mun track not available)";
0083 }
0084 return out;
0085 }
0086
0087 std::ostream& operator<<(std::ostream& out, TrackingParticleRef ref) {
0088 out << " [" << std::setw(4) << ref.key() << "]"
0089 << " type:" << std::setw(6) << ref->pdgId() << " pT: " << std::setw(6) << std::setprecision(3) << ref->pt()
0090 << " eta: " << std::setw(6) << std::setprecision(3) << ref->eta() << " phi: " << std::setw(6)
0091 << std::setprecision(3) << ref->phi();
0092 return out;
0093 }
0094
0095 void printAssociations(const char* label,
0096 TrackingParticleRef tp,
0097 const reco::SimToRecoCollection& byhits,
0098 const reco::SimToRecoCollection& bychi2) {
0099 reco::SimToRecoCollection::result_type found_byhits;
0100 if (byhits.find(tp) != byhits.end())
0101 found_byhits = byhits[tp];
0102 reco::SimToRecoCollection::result_type found_bychi2;
0103 if (bychi2.find(tp) != bychi2.end())
0104 found_bychi2 = bychi2[tp];
0105
0106 typedef std::tuple<double, double> Quality;
0107 std::map<edm::RefToBase<reco::Track>, Quality> found;
0108 for (std::vector<std::pair<edm::RefToBase<reco::Track>, double> >::const_iterator it = found_byhits.begin();
0109 it != found_byhits.end();
0110 ++it) {
0111 const edm::RefToBase<reco::Track> ref = it->first;
0112 found.insert(std::make_pair(ref, Quality()));
0113 std::get<0>(found[ref]) = it->second;
0114 }
0115 for (std::vector<std::pair<edm::RefToBase<reco::Track>, double> >::const_iterator it = found_bychi2.begin();
0116 it != found_bychi2.end();
0117 ++it) {
0118 const edm::RefToBase<reco::Track> ref = it->first;
0119 found.insert(std::make_pair(ref, Quality()));
0120 std::get<1>(found[ref]) = -it->second;
0121 }
0122
0123 for (std::map<edm::RefToBase<reco::Track>, Quality>::const_iterator it = found.begin(); it != found.end(); ++it) {
0124 std::cout << " " << std::setw(7) << std::left << label << std::right << it->first;
0125 if (std::get<0>(it->second))
0126 std::cout << " [" << std::setw(6) << std::setprecision(3) << std::get<0>(it->second) << "]";
0127 else
0128 std::cout << " ";
0129 if (std::get<1>(it->second))
0130 std::cout << " [" << std::setw(6) << std::setprecision(3) << std::get<1>(it->second) << "]";
0131 else
0132 std::cout << " ";
0133 std::cout << std::endl;
0134 }
0135 }
0136
0137 void printAssociations(const char* label,
0138 edm::RefToBase<reco::Track> tp,
0139 const reco::RecoToSimCollection& byhits,
0140 const reco::RecoToSimCollection& bychi2) {
0141 reco::RecoToSimCollection::result_type found_byhits;
0142 if (byhits.find(tp) != byhits.end())
0143 found_byhits = byhits[tp];
0144 reco::RecoToSimCollection::result_type found_bychi2;
0145 if (bychi2.find(tp) != bychi2.end())
0146 found_bychi2 = bychi2[tp];
0147
0148 typedef std::tuple<double, double> Quality;
0149 std::map<TrackingParticleRef, Quality> found;
0150 for (std::vector<std::pair<TrackingParticleRef, double> >::const_iterator it = found_byhits.begin();
0151 it != found_byhits.end();
0152 ++it) {
0153 const TrackingParticleRef ref = it->first;
0154 found.insert(std::make_pair(ref, Quality()));
0155 std::get<0>(found[ref]) = it->second;
0156 }
0157 for (std::vector<std::pair<TrackingParticleRef, double> >::const_iterator it = found_bychi2.begin();
0158 it != found_bychi2.end();
0159 ++it) {
0160 const TrackingParticleRef ref = it->first;
0161 found.insert(std::make_pair(ref, Quality()));
0162 std::get<1>(found[ref]) = -it->second;
0163 }
0164
0165 for (std::map<TrackingParticleRef, Quality>::const_iterator it = found.begin(); it != found.end(); ++it) {
0166 std::cout << " " << std::setw(7) << std::left << label << std::right << it->first;
0167 if (std::get<0>(it->second))
0168 std::cout << " [" << std::setw(6) << std::setprecision(3) << std::get<0>(it->second) << "]";
0169 else
0170 std::cout << " ";
0171 if (std::get<1>(it->second))
0172 std::cout << " [" << std::setw(6) << std::setprecision(3) << std::get<1>(it->second) << "]";
0173 else
0174 std::cout << " ";
0175 std::cout << std::endl;
0176 }
0177 }
0178
0179 void printAssociations(const char* label,
0180 reco::TrackRef tp,
0181 const reco::RecoToSimCollection& byhits,
0182 const reco::RecoToSimCollection& bychi2) {
0183 printAssociations(label, edm::RefToBase<reco::Track>(tp), byhits, bychi2);
0184 }
0185
0186 template <typename T>
0187 std::ostream& operator<<(std::ostream& out, const std::pair<edm::Ref<T>, double>& assoc) {
0188 out << assoc.first << " quality: " << assoc.second;
0189 return out;
0190 }
0191
0192 testLeptonAssociator::testLeptonAssociator(edm::ParameterSet const& iConfig) {
0193 token_recoTracks = consumes<edm::View<reco::Track> >(iConfig.getParameter<edm::InputTag>("tracks"));
0194 token_standAloneMuons =
0195 consumes<edm::View<reco::Track> >(iConfig.getParameter<edm::InputTag>("standAloneMuonTracks"));
0196 token_globalMuons = consumes<edm::View<reco::Track> >(iConfig.getParameter<edm::InputTag>("globalMuonTracks"));
0197 token_muons = consumes<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons"));
0198 token_trackingTruth = consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("trackingTruth"));
0199 m_flavour = iConfig.getParameter<unsigned int>("leptonFlavour");
0200 m_ptcut = iConfig.getParameter<double>("minPt");
0201 token_associatorByHits = consumes<reco::TrackToTrackingParticleAssociator>(edm::InputTag("trackAssociatorByHits"));
0202 token_associatorByChi2 = consumes<reco::TrackToTrackingParticleAssociator>(edm::InputTag("trackAssociatorByChi2"));
0203 }
0204
0205 void testLeptonAssociator::analyze(const edm::Event& iEvent, const edm::EventSetup& setup) {
0206 edm::Handle<edm::View<reco::Track> > recoTrackHandle;
0207 iEvent.getByToken(token_recoTracks, recoTrackHandle);
0208 const edm::View<reco::Track>& recoTrackCollection = *(recoTrackHandle.product());
0209
0210 edm::Handle<edm::View<reco::Track> > standAloneMuonHandle;
0211 iEvent.getByToken(token_standAloneMuons, standAloneMuonHandle);
0212 const edm::View<reco::Track>& standAloneMuonCollection = *(standAloneMuonHandle.product());
0213
0214 edm::Handle<edm::View<reco::Track> > globalMuonTrackHandle;
0215 iEvent.getByToken(token_globalMuons, globalMuonTrackHandle);
0216
0217
0218 edm::Handle<reco::MuonCollection> globalMuonHandle;
0219 iEvent.getByToken(token_muons, globalMuonHandle);
0220 const reco::MuonCollection& globalMuonCollection = *(globalMuonHandle.product());
0221
0222 edm::Handle<TrackingParticleCollection> trackingParticleHandle;
0223 iEvent.getByToken(token_trackingTruth, trackingParticleHandle);
0224 const TrackingParticleCollection& trackingParticleCollection = *(trackingParticleHandle.product());
0225
0226 edm::Handle<reco::TrackToTrackingParticleAssociator> associatorByHits;
0227 iEvent.getByToken(token_associatorByHits, associatorByHits);
0228
0229 edm::Handle<reco::TrackToTrackingParticleAssociator> associatorByChi2;
0230 iEvent.getByToken(token_associatorByChi2, associatorByChi2);
0231
0232 std::cout << std::fixed;
0233
0234 std::cout << std::endl;
0235 std::cout << "Found " << std::setw(6) << trackingParticleCollection.size() << " TrackingParticles" << std::flush;
0236 unsigned int count = 0;
0237 for (TrackingParticleCollection::size_type i = 0; i < trackingParticleCollection.size(); ++i)
0238 if ((std::abs(trackingParticleCollection[i].pdgId()) == (int)m_flavour) and
0239 (trackingParticleCollection[i].pt() >= m_ptcut))
0240 ++count;
0241
0242 std::cout << " ( " << std::setw(2) << count << " leptons with pT above " << m_ptcut << " GeV)" << std::endl;
0243 std::cout << " " << std::setw(6) << recoTrackCollection.size() << " Tracks" << std::endl;
0244 std::cout << " " << std::setw(6) << globalMuonCollection.size() << " Global muons" << std::endl;
0245 std::cout << " " << std::setw(6) << standAloneMuonCollection.size() << " StandAlone muons" << std::endl;
0246
0247
0248 {
0249 reco::SimToRecoCollection bychi2_tracks =
0250 associatorByChi2->associateSimToReco(recoTrackHandle, trackingParticleHandle);
0251 reco::SimToRecoCollection bychi2_globaltrack =
0252 associatorByChi2->associateSimToReco(globalMuonTrackHandle, trackingParticleHandle);
0253 reco::SimToRecoCollection bychi2_standalone =
0254 associatorByChi2->associateSimToReco(standAloneMuonHandle, trackingParticleHandle);
0255 reco::SimToRecoCollection byhits_tracks =
0256 associatorByHits->associateSimToReco(recoTrackHandle, trackingParticleHandle);
0257 reco::SimToRecoCollection byhits_globaltrack =
0258 associatorByHits->associateSimToReco(globalMuonTrackHandle, trackingParticleHandle);
0259 reco::SimToRecoCollection byhits_standalone =
0260 associatorByHits->associateSimToReco(standAloneMuonHandle, trackingParticleHandle);
0261
0262 for (TrackingParticleCollection::size_type i = 0; i < trackingParticleCollection.size(); ++i) {
0263 TrackingParticleRef tp(trackingParticleHandle, i);
0264 if ((std::abs(tp->pdgId()) != (int)m_flavour) or (tp->pt() < m_ptcut))
0265 continue;
0266 std::cout << "--> TrackingParticle" << tp << std::endl;
0267 printAssociations("Track", tp, byhits_tracks, bychi2_tracks);
0268 printAssociations("Local", tp, byhits_standalone, bychi2_standalone);
0269 printAssociations("Global", tp, byhits_globaltrack, bychi2_globaltrack);
0270 }
0271 }
0272
0273
0274 reco::RecoToSimCollection byhits_globalfake =
0275 associatorByHits->associateRecoToSim(recoTrackHandle, trackingParticleHandle);
0276 reco::RecoToSimCollection bychi2_globalfake =
0277 associatorByChi2->associateRecoToSim(recoTrackHandle, trackingParticleHandle);
0278 for (reco::MuonCollection::size_type i = 0; i < globalMuonCollection.size(); ++i) {
0279 reco::MuonRef lepton(globalMuonHandle, i);
0280 std::cout << "<-- Global " << lepton << std::endl;
0281 printAssociations("TrackingParticle", lepton->track(), byhits_globalfake, bychi2_globalfake);
0282 }
0283
0284
0285 reco::RecoToSimCollection byhits_standalonefake =
0286 associatorByHits->associateRecoToSim(standAloneMuonHandle, trackingParticleHandle);
0287 reco::RecoToSimCollection bychi2_standalonefake =
0288 associatorByChi2->associateRecoToSim(standAloneMuonHandle, trackingParticleHandle);
0289 for (edm::View<reco::Track>::size_type i = 0; i < standAloneMuonCollection.size(); ++i) {
0290 edm::RefToBase<reco::Track> lepton(standAloneMuonHandle, i);
0291 std::cout << "<-- Local " << lepton << std::endl;
0292 printAssociations("TrackingParticle", lepton, byhits_standalonefake, bychi2_standalonefake);
0293 }
0294 }
0295
0296 #include "FWCore/Framework/interface/MakerMacros.h"
0297 DEFINE_FWK_MODULE(testLeptonAssociator);