RecoMuon/src/GlobalMuonMatchAnalyzer.cc

0001 /**
0002  * Class: GlobalMuonMatchAnalyzer
0003  *
0004  *
0005  *
0006  * Authors :
0007  * \author Adam Everett - Purdue University
0008  *
0009  */
0010
0011 #include "Validation/RecoMuon/src/GlobalMuonMatchAnalyzer.h"
0012
0013 // user include files
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021
0022 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0023
0024 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0025 #include "DataFormats/MuonReco/interface/MuonTrackLinks.h"
0026 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0027
0028 #include "DQMServices/Core/interface/DQMStore.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030
0031 #include <TH2.h>
0032
0033 GlobalMuonMatchAnalyzer::GlobalMuonMatchAnalyzer(const edm::ParameterSet &ps)
0034
0035 {
0036   iConfig = ps;
0037   //now do what ever initialization is needed
0038   tkAssociatorName_ = iConfig.getUntrackedParameter<edm::InputTag>("tkAssociator");
0039   muAssociatorName_ = iConfig.getUntrackedParameter<edm::InputTag>("muAssociator");
0040
0041   tkAssociatorToken_ = consumes<reco::TrackToTrackingParticleAssociator>(tkAssociatorName_);
0042   muAssociatorToken_ = consumes<reco::TrackToTrackingParticleAssociator>(muAssociatorName_);
0043
0044   subsystemname_ = iConfig.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
0045   tpName_ = iConfig.getUntrackedParameter<edm::InputTag>("tpLabel");
0046   tkName_ = iConfig.getUntrackedParameter<edm::InputTag>("tkLabel");
0047   staName_ = iConfig.getUntrackedParameter<edm::InputTag>("muLabel");
0048   glbName_ = iConfig.getUntrackedParameter<edm::InputTag>("glbLabel");
0049
0050   out = iConfig.getUntrackedParameter<std::string>("out");
0051   dbe_ = edm::Service<DQMStore>().operator->();
0052
0053   tpToken_ = consumes<edm::View<reco::Track> >(tpName_);
0054   tkToken_ = consumes<edm::View<reco::Track> >(tkName_);
0055   staToken_ = consumes<edm::View<reco::Track> >(staName_);
0056   glbToken_ = consumes<edm::View<reco::Track> >(glbName_);
0057 }
0058
0059 GlobalMuonMatchAnalyzer::~GlobalMuonMatchAnalyzer() {
0060   // do anything here that needs to be done at desctruction time
0061   // (e.g. close files, deallocate resources etc.)
0062 }
0063
0064 //
0065 // member functions
0066 //
0067
0068 // ------------ method called to for each event  ------------
0069 void GlobalMuonMatchAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0070   using namespace edm;
0071   using namespace reco;
0072
0073   Handle<TrackingParticleCollection> tpHandle;
0074   iEvent.getByToken(tpToken_, tpHandle);
0075   const TrackingParticleCollection tpColl = *(tpHandle.product());
0076
0077   Handle<reco::MuonTrackLinksCollection> muHandle;
0078   iEvent.getByToken(glbToken_, muHandle);
0079   const reco::MuonTrackLinksCollection muColl = *(muHandle.product());
0080
0081   Handle<View<Track> > staHandle;
0082   iEvent.getByToken(staToken_, staHandle);
0083
0084   Handle<View<Track> > glbHandle;
0085   iEvent.getByToken(glbToken_, glbHandle);
0086
0087   Handle<View<Track> > tkHandle;
0088   iEvent.getByToken(tkToken_, tkHandle);
0089
0090   edm::Handle<reco::TrackToTrackingParticleAssociator> tkAssociator;
0091   iEvent.getByToken(tkAssociatorToken_, tkAssociator);
0092
0093   // Mu Associator
0094   edm::Handle<reco::TrackToTrackingParticleAssociator> muAssociator;
0095   iEvent.getByToken(muAssociatorToken_, muAssociator);
0096
0097   reco::RecoToSimCollection tkrecoToSimCollection = tkAssociator->associateRecoToSim(tkHandle, tpHandle);
0098   reco::SimToRecoCollection tksimToRecoCollection = tkAssociator->associateSimToReco(tkHandle, tpHandle);
0099
0100   reco::RecoToSimCollection starecoToSimCollection = muAssociator->associateRecoToSim(staHandle, tpHandle);
0101   reco::SimToRecoCollection stasimToRecoCollection = muAssociator->associateSimToReco(staHandle, tpHandle);
0102
0103   reco::RecoToSimCollection glbrecoToSimCollection = muAssociator->associateRecoToSim(glbHandle, tpHandle);
0104   reco::SimToRecoCollection glbsimToRecoCollection = muAssociator->associateSimToReco(glbHandle, tpHandle);
0105
0106   for (TrackingParticleCollection::size_type i = 0; i < tpColl.size(); ++i) {
0107     TrackingParticleRef tp(tpHandle, i);
0108
0109     std::vector<std::pair<RefToBase<Track>, double> > rvGlb;
0110     RefToBase<Track> rGlb;
0111     if (glbsimToRecoCollection.find(tp) != glbsimToRecoCollection.end()) {
0112       rvGlb = glbsimToRecoCollection[tp];
0113       if (!rvGlb.empty()) {
0114         rGlb = rvGlb.begin()->first;
0115       }
0116     }
0117
0118     std::vector<std::pair<RefToBase<Track>, double> > rvSta;
0119     RefToBase<Track> rSta;
0120     if (stasimToRecoCollection.find(tp) != stasimToRecoCollection.end()) {
0121       rvSta = stasimToRecoCollection[tp];
0122       if (!rvSta.empty()) {
0123         rSta = rvSta.begin()->first;
0124       }
0125     }
0126
0127     std::vector<std::pair<RefToBase<Track>, double> > rvTk;
0128     RefToBase<Track> rTk;
0129     if (tksimToRecoCollection.find(tp) != tksimToRecoCollection.end()) {
0130       rvTk = tksimToRecoCollection[tp];
0131       if (!rvTk.empty()) {
0132         rTk = rvTk.begin()->first;
0133       }
0134     }
0135
0136     if (!rvSta.empty() && !rvTk.empty()) {
0137       //should have matched
0138       h_shouldMatch->Fill(rTk->eta(), rTk->pt());
0139     }
0140
0141     for (reco::MuonTrackLinksCollection::const_iterator links = muHandle->begin(); links != muHandle->end(); ++links) {
0142       if (rGlb == RefToBase<Track>(links->globalTrack())) {
0143         if (RefToBase<Track>(links->trackerTrack()) == rTk && RefToBase<Track>(links->standAloneTrack()) == rSta) {
0144           //goodMatch
0145           h_goodMatchSim->Fill(rGlb->eta(), rGlb->pt());
0146         }
0147         if (RefToBase<Track>(links->trackerTrack()) == rTk && RefToBase<Track>(links->standAloneTrack()) != rSta) {
0148           //tkOnlyMatch
0149           h_tkOnlySim->Fill(rGlb->eta(), rGlb->pt());
0150         }
0151         if (RefToBase<Track>(links->standAloneTrack()) == rSta && RefToBase<Track>(links->trackerTrack()) != rTk) {
0152           //staOnlyMatch
0153           h_staOnlySim->Fill(rGlb->eta(), rGlb->pt());
0154         }
0155       }
0156     }
0157   }
0158
0159   ////////
0160
0161   for (reco::MuonTrackLinksCollection::const_iterator links = muHandle->begin(); links != muHandle->end(); ++links) {
0162     RefToBase<Track> glbRef = RefToBase<Track>(links->globalTrack());
0163     RefToBase<Track> staRef = RefToBase<Track>(links->standAloneTrack());
0164     RefToBase<Track> tkRef = RefToBase<Track>(links->trackerTrack());
0165
0166     std::vector<std::pair<TrackingParticleRef, double> > tp1;
0167     TrackingParticleRef tp1r;
0168     if (glbrecoToSimCollection.find(glbRef) != glbrecoToSimCollection.end()) {
0169       tp1 = glbrecoToSimCollection[glbRef];
0170       if (!tp1.empty()) {
0171         tp1r = tp1.begin()->first;
0172       }
0173     }
0174
0175     std::vector<std::pair<TrackingParticleRef, double> > tp2;
0176     TrackingParticleRef tp2r;
0177     if (starecoToSimCollection.find(staRef) != starecoToSimCollection.end()) {
0178       tp2 = starecoToSimCollection[staRef];
0179       if (!tp2.empty()) {
0180         tp2r = tp2.begin()->first;
0181       }
0182     }
0183
0184     std::vector<std::pair<TrackingParticleRef, double> > tp3;
0185     TrackingParticleRef tp3r;
0186     if (tkrecoToSimCollection.find(tkRef) != tkrecoToSimCollection.end()) {
0187       tp3 = tkrecoToSimCollection[tkRef];
0188       if (!tp3.empty()) {
0189         tp3r = tp3.begin()->first;
0190       }
0191     }
0192
0193     if (!tp1.empty()) {
0194       //was reconstructed
0195       h_totReco->Fill(glbRef->eta(), glbRef->pt());
0196       if (tp2r == tp3r) {  // && tp1r == tp3r) {
0197         //came from same TP
0198         h_goodMatch->Fill(glbRef->eta(), glbRef->pt());
0199       } else {
0200         //mis-match
0201         h_fakeMatch->Fill(glbRef->eta(), glbRef->pt());
0202       }
0203     }
0204   }
0205 }
0206
0207 // ------------ method called once each job just after ending the event loop  ------------
0208 void GlobalMuonMatchAnalyzer::dqmEndRun(edm::Run const &, edm::EventSetup const &) {
0209   computeEfficiencyEta(h_effic, h_goodMatchSim, h_shouldMatch);
0210   computeEfficiencyPt(h_efficPt, h_goodMatchSim, h_shouldMatch);
0211
0212   computeEfficiencyEta(h_fake, h_fakeMatch, h_totReco);
0213   computeEfficiencyPt(h_fakePt, h_fakeMatch, h_totReco);
0214
0215   if (!out.empty() && dbe_)
0216     dbe_->save(out);
0217 }
0218
0219 //void GlobalMuonMatchAnalyzer::beginRun(const edm::Run&, const edm::EventSetup& setup)
0220
0221 void GlobalMuonMatchAnalyzer::bookHistograms(DQMStore::IBooker &ibooker,
0222                                              edm::Run const &iRun,
0223                                              edm::EventSetup const &iSetup) {
0224   // Tk Associator
0225
0226   ibooker.cd();
0227   // Run histos only for dqmEndRun handling
0228   ibooker.setScope(MonitorElementData::Scope::RUN);
0229
0230   std::string dirName = "Matcher/";
0231   //  ibooker.setCurrentFolder("RecoMuonV/Matcher");
0232   ibooker.setCurrentFolder(dirName);
0233
0234   h_shouldMatch = ibooker.book2D("h_shouldMatch", "SIM associated to Tk and Sta", 50, -2.5, 2.5, 100, 0., 500.);
0235   h_goodMatchSim = ibooker.book2D("h_goodMatchSim", "SIM associated to Glb Sta Tk", 50, -2.5, 2.5, 100, 0., 500.);
0236   h_tkOnlySim = ibooker.book2D("h_tkOnlySim", "SIM associated to Glb Tk", 50, -2.5, 2.5, 100, 0., 500.);
0237   h_staOnlySim = ibooker.book2D("h_staOnlySim", "SIM associated to Glb Sta", 50, -2.5, 2.5, 100, 0., 500.);
0238
0239   h_totReco = ibooker.book2D("h_totReco", "Total Glb Reconstructed", 50, -2.5, 2.5, 100, 0., 500.);
0240   h_goodMatch = ibooker.book2D("h_goodMatch", "Sta and Tk from same SIM", 50, -2.5, 2.5, 100, 0., 500.);
0241   h_fakeMatch = ibooker.book2D("h_fakeMatch", "Sta and Tk not from same SIM", 50, -2.5, 2.5, 100, 0., 500.);
0242
0243   h_effic = ibooker.book1D("h_effic", "Efficiency vs #eta", 50, -2.5, 2.5);
0244   h_efficPt = ibooker.book1D("h_efficPt", "Efficiency vs p_{T}", 100, 0., 100.);
0245
0246   h_fake = ibooker.book1D("h_fake", "Fake fraction vs #eta", 50, -2.5, 2.5);
0247   h_fakePt = ibooker.book1D("h_fakePt", "Fake fraction vs p_{T}", 100, 0., 100.);
0248 }
0249
0250 void GlobalMuonMatchAnalyzer::computeEfficiencyEta(MonitorElement *effHist,
0251                                                    MonitorElement *recoTH2,
0252                                                    MonitorElement *simTH2) {
0253   TH2F *h1 = recoTH2->getTH2F();
0254   TH1D *reco = h1->ProjectionX();
0255
0256   TH2F *h2 = simTH2->getTH2F();
0257   TH1D *sim = h2->ProjectionX();
0258
0259   TH1F *hEff = (TH1F *)reco->Clone();
0260
0261   hEff->Divide(sim);
0262
0263   hEff->SetName("tmp_" + TString(reco->GetName()));
0264
0265   // Set the error accordingly to binomial statistics
0266   int nBinsEta = hEff->GetNbinsX();
0267   for (int bin = 1; bin <= nBinsEta; bin++) {
0268     float nSimHit = sim->GetBinContent(bin);
0269     float eff = hEff->GetBinContent(bin);
0270     float error = 0;
0271     if (nSimHit != 0 && eff <= 1) {
0272       error = sqrt(eff * (1 - eff) / nSimHit);
0273     }
0274     hEff->SetBinError(bin, error);
0275     effHist->setBinContent(bin, eff);
0276     effHist->setBinError(bin, error);
0277   }
0278 }
0279
0280 void GlobalMuonMatchAnalyzer::computeEfficiencyPt(MonitorElement *effHist,
0281                                                   MonitorElement *recoTH2,
0282                                                   MonitorElement *simTH2) {
0283   TH2F *h1 = recoTH2->getTH2F();
0284   TH1D *reco = h1->ProjectionY();
0285
0286   TH2F *h2 = simTH2->getTH2F();
0287   TH1D *sim = h2->ProjectionY();
0288
0289   TH1F *hEff = (TH1F *)reco->Clone();
0290
0291   hEff->Divide(sim);
0292
0293   hEff->SetName("tmp_" + TString(reco->GetName()));
0294
0295   // Set the error accordingly to binomial statistics
0296   int nBinsPt = hEff->GetNbinsX();
0297   for (int bin = 1; bin <= nBinsPt; bin++) {
0298     float nSimHit = sim->GetBinContent(bin);
0299     float eff = hEff->GetBinContent(bin);
0300     float error = 0;
0301     if (nSimHit != 0 && eff <= 1) {
0302       error = sqrt(eff * (1 - eff) / nSimHit);
0303     }
0304     hEff->SetBinError(bin, error);
0305     effHist->setBinContent(bin, eff);
0306     effHist->setBinError(bin, error);
0307   }
0308 }