File indexing completed on 2023-03-17 11:28:34
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "Validation/RecoMuon/src/GlobalMuonMatchAnalyzer.h"
0012
0013
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
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
0061
0062 }
0063
0064
0065
0066
0067
0068
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
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
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
0145 h_goodMatchSim->Fill(rGlb->eta(), rGlb->pt());
0146 }
0147 if (RefToBase<Track>(links->trackerTrack()) == rTk && RefToBase<Track>(links->standAloneTrack()) != rSta) {
0148
0149 h_tkOnlySim->Fill(rGlb->eta(), rGlb->pt());
0150 }
0151 if (RefToBase<Track>(links->standAloneTrack()) == rSta && RefToBase<Track>(links->trackerTrack()) != rTk) {
0152
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
0195 h_totReco->Fill(glbRef->eta(), glbRef->pt());
0196 if (tp2r == tp3r) {
0197
0198 h_goodMatch->Fill(glbRef->eta(), glbRef->pt());
0199 } else {
0200
0201 h_fakeMatch->Fill(glbRef->eta(), glbRef->pt());
0202 }
0203 }
0204 }
0205 }
0206
0207
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
0220
0221 void GlobalMuonMatchAnalyzer::bookHistograms(DQMStore::IBooker &ibooker,
0222 edm::Run const &iRun,
0223 edm::EventSetup const &iSetup) {
0224
0225
0226 ibooker.cd();
0227
0228 ibooker.setScope(MonitorElementData::Scope::RUN);
0229
0230 std::string dirName = "Matcher/";
0231
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
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
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 }