Back to home page

Project CMSSW displayed by LXR
 
 

    

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 }  // namespace std
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;  // why is chi2 negative ?
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;  // why is chi2 negative ?
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   //const edm::View<reco::Track>& globalMuonTrackCollection = *(globalMuonTrackHandle.product());
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   // look for tracks and muons associated to the tracking particles
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   // look for tracking particles associated to the (tracker part of the) reconstructed global muons
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   // look for tracking particles associated to the reconstructed standAlone muons
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);
[ 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!