Back to home page

Project CMSSW displayed by LXR
 
 

    

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 }  // namespace std
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   // access EventSetup during the Event loop, to make sure conditions are up to date
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   // access Event collections
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   // look for tracks and muons associated to the tracking particles
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     try {
0294       bypos_tracks       = associatorByPos ->associateSimToReco(recoTrackHandle,       trackingParticleHandle );
0295     } catch (const std::exception & e) { std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl; }
0296     try {
0297       bypos_globaltrack  = associatorByPos ->associateSimToReco(globalMuonTrackHandle, trackingParticleHandle );
0298     } catch (const std::exception & e) { std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl; }
0299     try {
0300       bypos_standalone   = associatorByPos ->associateSimToReco(standAloneMuonHandle,  trackingParticleHandle );
0301     } catch (const std::exception & e) { std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl; }
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   // look for tracking particles associated to the (tracker part of the) reconstructed global muons
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   try {
0333     bypos_muon  = associatorByPos ->associateRecoToSim (recoTrackHandle, trackingParticleHandle );
0334   } catch (const std::exception & e) { std::cerr << std::endl << "Error building reco::SimToRecoCollection:\n" << e.what() << std::endl; }
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   // look for tracking particles associated to the reconstructed Global muons
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   try {
0359     bypos_global  = associatorByPos ->associateRecoToSim (globalMuonTrackHandle, trackingParticleHandle );
0360   } catch (const std::exception & e) { std::cerr << std::endl << "Error building reco::RecoToSimCollection:\n" << e.what() << std::endl; }
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   // look for tracking particles associated to the reconstructed StandAlone muons
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   try {
0385     bypos_standalone  = associatorByPos ->associateRecoToSim (standAloneMuonHandle, trackingParticleHandle );
0386   } catch (const std::exception & e) { std::cerr << std::endl << "Error building reco::RecoToSimCollection:\n" << e.what() << std::endl; }
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!