File indexing completed on 2024-08-13 22:17:12
0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "Validation/RPCRecHits/interface/RPCRecHitValid.h"
0003
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0008 #include "DataFormats/MuonReco/interface/Muon.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0011 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0012 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0013 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
0014 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0015 #include "Geometry/RPCGeometry/interface/RPCRollSpecs.h"
0016 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0017
0018 #include <algorithm>
0019
0020 using namespace std;
0021
0022 typedef RPCRecHitValid::MonitorElement *MEP;
0023
0024 RPCRecHitValid::RPCRecHitValid(const edm::ParameterSet &pset) {
0025 simHitToken_ = consumes<SimHits>(pset.getParameter<edm::InputTag>("simHit"));
0026 recHitToken_ = consumes<RecHits>(pset.getParameter<edm::InputTag>("recHit"));
0027 simParticleToken_ = consumes<SimParticles>(pset.getParameter<edm::InputTag>("simTrack"));
0028 simHitAssocToken_ = consumes<SimHitAssoc>(pset.getParameter<edm::InputTag>("simHitAssoc"));
0029 muonToken_ = consumes<reco::MuonCollection>(pset.getParameter<edm::InputTag>("muon"));
0030
0031 subDir_ = pset.getParameter<std::string>("subDir");
0032
0033 rpcGeomToken_ = esConsumes();
0034 rpcGeomTokenInRun_ = esConsumes<edm::Transition::BeginRun>();
0035 }
0036
0037 void RPCRecHitValid::bookHistograms(DQMStore::IBooker &booker, edm::Run const &run, edm::EventSetup const &eventSetup) {
0038
0039 h_.bookHistograms(booker, subDir_);
0040
0041
0042 booker.setCurrentFolder(subDir_ + "/HitProperty");
0043 h_simParticleType = booker.book1D("SimHitPType", "SimHit particle type", 11, 0, 11);
0044 h_simParticleType->getTH1()->SetMinimum(0);
0045 if (TH1 *h = h_simParticleType->getTH1()) {
0046 h->GetXaxis()->SetBinLabel(1, "#mu^{-}");
0047 h->GetXaxis()->SetBinLabel(2, "#mu^{+}");
0048 h->GetXaxis()->SetBinLabel(3, "e^{-}");
0049 h->GetXaxis()->SetBinLabel(4, "e^{+}");
0050 h->GetXaxis()->SetBinLabel(5, "#pi^{+}");
0051 h->GetXaxis()->SetBinLabel(6, "#pi^{-}");
0052 h->GetXaxis()->SetBinLabel(7, "K^{+}");
0053 h->GetXaxis()->SetBinLabel(8, "K^{-}");
0054 h->GetXaxis()->SetBinLabel(9, "p^{+}");
0055 h->GetXaxis()->SetBinLabel(10, "p^{-}");
0056 h->GetXaxis()->SetBinLabel(11, "Other");
0057 }
0058
0059 booker.setCurrentFolder(subDir_ + "/Track");
0060
0061 h_nRPCHitPerSimMuon = booker.book1D("NRPCHitPerSimMuon", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
0062 h_nRPCHitPerSimMuonBarrel =
0063 booker.book1D("NRPCHitPerSimMuonBarrel", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
0064 h_nRPCHitPerSimMuonOverlap =
0065 booker.book1D("NRPCHitPerSimMuonOverlap", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
0066 h_nRPCHitPerSimMuonEndcap =
0067 booker.book1D("NRPCHitPerSimMuonEndcap", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
0068
0069 h_nRPCHitPerRecoMuon = booker.book1D("NRPCHitPerRecoMuon", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
0070 h_nRPCHitPerRecoMuonBarrel =
0071 booker.book1D("NRPCHitPerRecoMuonBarrel", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
0072 h_nRPCHitPerRecoMuonOverlap =
0073 booker.book1D("NRPCHitPerRecoMuonOverlap", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
0074 h_nRPCHitPerRecoMuonEndcap =
0075 booker.book1D("NRPCHitPerRecoMuonEndcap", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
0076
0077 h_nRPCHitPerSimMuon->getTH1()->SetMinimum(0);
0078 h_nRPCHitPerSimMuonBarrel->getTH1()->SetMinimum(0);
0079 h_nRPCHitPerSimMuonOverlap->getTH1()->SetMinimum(0);
0080 h_nRPCHitPerSimMuonEndcap->getTH1()->SetMinimum(0);
0081
0082 h_nRPCHitPerRecoMuon->getTH1()->SetMinimum(0);
0083 h_nRPCHitPerRecoMuonBarrel->getTH1()->SetMinimum(0);
0084 h_nRPCHitPerRecoMuonOverlap->getTH1()->SetMinimum(0);
0085 h_nRPCHitPerRecoMuonEndcap->getTH1()->SetMinimum(0);
0086
0087 float ptBins[] = {0, 1, 2, 5, 10, 20, 30, 50, 100, 200, 300, 500};
0088 const int nPtBins = sizeof(ptBins) / sizeof(float) - 1;
0089 h_simMuonBarrel_pt =
0090 booker.book1D("SimMuonBarrel_pt", "SimMuon RPCHit in Barrel p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0091 h_simMuonOverlap_pt =
0092 booker.book1D("SimMuonOverlap_pt", "SimMuon RPCHit in Overlap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0093 h_simMuonEndcap_pt =
0094 booker.book1D("SimMuonEndcap_pt", "SimMuon RPCHit in Endcap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0095 h_simMuonNoRPC_pt =
0096 booker.book1D("SimMuonNoRPC_pt", "SimMuon without RPCHit p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0097 h_simMuonBarrel_eta = booker.book1D("SimMuonBarrel_eta", "SimMuon RPCHit in Barrel #eta;#eta", 50, -2.5, 2.5);
0098 h_simMuonOverlap_eta = booker.book1D("SimMuonOverlap_eta", "SimMuon RPCHit in Overlap #eta;#eta", 50, -2.5, 2.5);
0099 h_simMuonEndcap_eta = booker.book1D("SimMuonEndcap_eta", "SimMuon RPCHit in Endcap #eta;#eta", 50, -2.5, 2.5);
0100 h_simMuonNoRPC_eta = booker.book1D("SimMuonNoRPC_eta", "SimMuon without RPCHit #eta;#eta", 50, -2.5, 2.5);
0101 h_simMuonBarrel_phi =
0102 booker.book1D("SimMuonBarrel_phi", "SimMuon RPCHit in Barrel #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0103 h_simMuonOverlap_phi =
0104 booker.book1D("SimMuonOverlap_phi", "SimMuon RPCHit in Overlap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0105 h_simMuonEndcap_phi =
0106 booker.book1D("SimMuonEndcap_phi", "SimMuon RPCHit in Endcap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0107 h_simMuonNoRPC_phi =
0108 booker.book1D("SimMuonNoRPC_phi", "SimMuon without RPCHit #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0109
0110 h_recoMuonBarrel_pt =
0111 booker.book1D("RecoMuonBarrel_pt", "RecoMuon RPCHit in Barrel p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0112 h_recoMuonOverlap_pt =
0113 booker.book1D("RecoMuonOverlap_pt", "RecoMuon RPCHit in Overlap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0114 h_recoMuonEndcap_pt =
0115 booker.book1D("RecoMuonEndcap_pt", "RecoMuon RPCHit in Endcap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0116 h_recoMuonNoRPC_pt =
0117 booker.book1D("RecoMuonNoRPC_pt", "RecoMuon without RPCHit p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
0118 h_recoMuonBarrel_eta = booker.book1D("RecoMuonBarrel_eta", "RecoMuon RPCHit in Barrel #eta;#eta", 50, -2.5, 2.5);
0119 h_recoMuonOverlap_eta = booker.book1D("RecoMuonOverlap_eta", "RecoMuon RPCHit in Overlap #eta;#eta", 50, -2.5, 2.5);
0120 h_recoMuonEndcap_eta = booker.book1D("RecoMuonEndcap_eta", "RecoMuon RPCHit in Endcap #eta;#eta", 50, -2.5, 2.5);
0121 h_recoMuonNoRPC_eta = booker.book1D("RecoMuonNoRPC_eta", "RecoMuon without RPCHit #eta;#eta", 50, -2.5, 2.5);
0122 h_recoMuonBarrel_phi =
0123 booker.book1D("RecoMuonBarrel_phi", "RecoMuon RPCHit in Barrel #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0124 h_recoMuonOverlap_phi =
0125 booker.book1D("RecoMuonOverlap_phi", "RecoMuon RPCHit in Overlap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0126 h_recoMuonEndcap_phi =
0127 booker.book1D("RecoMuonEndcap_phi", "RecoMuon RPCHit in Endcap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0128 h_recoMuonNoRPC_phi =
0129 booker.book1D("RecoMuonNoRPC_phi", "RecoMuon without RPCHit #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
0130
0131 h_simMuonBarrel_pt->getTH1()->SetMinimum(0);
0132 h_simMuonOverlap_pt->getTH1()->SetMinimum(0);
0133 h_simMuonEndcap_pt->getTH1()->SetMinimum(0);
0134 h_simMuonNoRPC_pt->getTH1()->SetMinimum(0);
0135 h_simMuonBarrel_eta->getTH1()->SetMinimum(0);
0136 h_simMuonOverlap_eta->getTH1()->SetMinimum(0);
0137 h_simMuonEndcap_eta->getTH1()->SetMinimum(0);
0138 h_simMuonNoRPC_eta->getTH1()->SetMinimum(0);
0139 h_simMuonBarrel_phi->getTH1()->SetMinimum(0);
0140 h_simMuonOverlap_phi->getTH1()->SetMinimum(0);
0141 h_simMuonEndcap_phi->getTH1()->SetMinimum(0);
0142 h_simMuonNoRPC_phi->getTH1()->SetMinimum(0);
0143
0144 h_recoMuonBarrel_pt->getTH1()->SetMinimum(0);
0145 h_recoMuonOverlap_pt->getTH1()->SetMinimum(0);
0146 h_recoMuonEndcap_pt->getTH1()->SetMinimum(0);
0147 h_recoMuonNoRPC_pt->getTH1()->SetMinimum(0);
0148 h_recoMuonBarrel_eta->getTH1()->SetMinimum(0);
0149 h_recoMuonOverlap_eta->getTH1()->SetMinimum(0);
0150 h_recoMuonEndcap_eta->getTH1()->SetMinimum(0);
0151 h_recoMuonNoRPC_eta->getTH1()->SetMinimum(0);
0152 h_recoMuonBarrel_phi->getTH1()->SetMinimum(0);
0153 h_recoMuonOverlap_phi->getTH1()->SetMinimum(0);
0154 h_recoMuonEndcap_phi->getTH1()->SetMinimum(0);
0155 h_recoMuonNoRPC_phi->getTH1()->SetMinimum(0);
0156
0157 booker.setCurrentFolder(subDir_ + "/Occupancy");
0158
0159 h_eventCount = booker.book1D("EventCount", "Event count", 3, 1, 4);
0160 h_eventCount->getTH1()->SetMinimum(0);
0161 if (h_eventCount) {
0162 TH1 *h = h_eventCount->getTH1();
0163 h->GetXaxis()->SetBinLabel(1, "eventBegin");
0164 h->GetXaxis()->SetBinLabel(2, "eventEnd");
0165 h->GetXaxis()->SetBinLabel(3, "run");
0166 }
0167 h_eventCount->Fill(3);
0168
0169
0170 auto rpcGeom = eventSetup.getHandle(rpcGeomTokenInRun_);
0171
0172 int nRPCRollBarrel = 0, nRPCRollEndcap = 0;
0173
0174 TrackingGeometry::DetContainer rpcDets = rpcGeom->dets();
0175 for (auto det : rpcDets) {
0176 auto rpcCh = dynamic_cast<const RPCChamber *>(det);
0177 if (!rpcCh)
0178 continue;
0179
0180 std::vector<const RPCRoll *> rolls = rpcCh->rolls();
0181 for (auto roll : rolls) {
0182 if (!roll)
0183 continue;
0184
0185 const int rawId = roll->geographicalId().rawId();
0186
0187 if (roll->isBarrel()) {
0188 detIdToIndexMapBarrel_[rawId] = nRPCRollBarrel;
0189 ++nRPCRollBarrel;
0190 } else {
0191 detIdToIndexMapEndcap_[rawId] = nRPCRollEndcap;
0192 ++nRPCRollEndcap;
0193 }
0194 }
0195 }
0196
0197 booker.setCurrentFolder(subDir_ + "/Occupancy");
0198 h_matchOccupancyBarrel_detId = booker.book1D("MatchOccupancyBarrel_detId",
0199 "Matched hit occupancy;roll index (can be arbitrary)",
0200 nRPCRollBarrel,
0201 0,
0202 nRPCRollBarrel);
0203 h_matchOccupancyEndcap_detId = booker.book1D("MatchOccupancyEndcap_detId",
0204 "Matched hit occupancy;roll index (can be arbitrary)",
0205 nRPCRollEndcap,
0206 0,
0207 nRPCRollEndcap);
0208 h_refOccupancyBarrel_detId = booker.book1D("RefOccupancyBarrel_detId",
0209 "Reference hit occupancy;roll index (can be arbitrary)",
0210 nRPCRollBarrel,
0211 0,
0212 nRPCRollBarrel);
0213 h_refOccupancyEndcap_detId = booker.book1D("RefOccupancyEndcap_detId",
0214 "Reference hit occupancy;roll index (can be arbitrary)",
0215 nRPCRollEndcap,
0216 0,
0217 nRPCRollEndcap);
0218 h_allOccupancyBarrel_detId = booker.book1D(
0219 "OccupancyBarrel_detId", "Occupancy;roll index (can be arbitrary)", nRPCRollBarrel, 0, nRPCRollBarrel);
0220 h_allOccupancyEndcap_detId = booker.book1D(
0221 "OccupancyEndcap_detId", "Occupancy;roll index (can be arbitrary)", nRPCRollEndcap, 0, nRPCRollEndcap);
0222
0223 h_matchOccupancyBarrel_detId->getTH1()->SetMinimum(0);
0224 h_matchOccupancyEndcap_detId->getTH1()->SetMinimum(0);
0225 h_refOccupancyBarrel_detId->getTH1()->SetMinimum(0);
0226 h_refOccupancyEndcap_detId->getTH1()->SetMinimum(0);
0227 h_allOccupancyBarrel_detId->getTH1()->SetMinimum(0);
0228 h_allOccupancyEndcap_detId->getTH1()->SetMinimum(0);
0229
0230 h_rollAreaBarrel_detId = booker.bookProfile(
0231 "RollAreaBarrel_detId", "Roll area;roll index;Area", nRPCRollBarrel, 0., 1. * nRPCRollBarrel, 0., 1e5);
0232 h_rollAreaEndcap_detId = booker.bookProfile(
0233 "RollAreaEndcap_detId", "Roll area;roll index;Area", nRPCRollEndcap, 0., 1. * nRPCRollEndcap, 0., 1e5);
0234
0235 for (auto detIdToIndex : detIdToIndexMapBarrel_) {
0236 const int rawId = detIdToIndex.first;
0237 const int index = detIdToIndex.second;
0238
0239 const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
0240 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
0241
0242 const StripTopology &topol = roll->specificTopology();
0243 const double area = topol.stripLength() * topol.nstrips() * topol.pitch();
0244
0245 h_rollAreaBarrel_detId->Fill(index, area);
0246 }
0247
0248 for (auto detIdToIndex : detIdToIndexMapEndcap_) {
0249 const int rawId = detIdToIndex.first;
0250 const int index = detIdToIndex.second;
0251
0252 const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
0253 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
0254
0255 const StripTopology &topol = roll->specificTopology();
0256 const double area = topol.stripLength() * topol.nstrips() * topol.pitch();
0257
0258 h_rollAreaEndcap_detId->Fill(index, area);
0259 }
0260 }
0261
0262 void RPCRecHitValid::analyze(const edm::Event &event, const edm::EventSetup &eventSetup) {
0263 h_eventCount->Fill(1);
0264
0265
0266 auto rpcGeom = eventSetup.getHandle(rpcGeomToken_);
0267
0268
0269 edm::Handle<edm::PSimHitContainer> simHitHandle;
0270 if (!event.getByToken(simHitToken_, simHitHandle)) {
0271 edm::LogInfo("RPCRecHitValid") << "Cannot find simHit collection\n";
0272 return;
0273 }
0274
0275
0276 edm::Handle<RPCRecHitCollection> recHitHandle;
0277 if (!event.getByToken(recHitToken_, recHitHandle)) {
0278 edm::LogInfo("RPCRecHitValid") << "Cannot find recHit collection\n";
0279 return;
0280 }
0281
0282
0283 edm::Handle<TrackingParticleCollection> simParticleHandle;
0284 if (!event.getByToken(simParticleToken_, simParticleHandle)) {
0285 edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle collection\n";
0286 return;
0287 }
0288
0289
0290 edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0291 if (!event.getByToken(simHitAssocToken_, simHitsTPAssoc)) {
0292 edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle to SimHit association map\n";
0293 return;
0294 }
0295
0296
0297 edm::Handle<reco::MuonCollection> muonHandle;
0298 if (!event.getByToken(muonToken_, muonHandle)) {
0299 edm::LogInfo("RPCRecHitValid") << "Cannot find muon collection\n";
0300 return;
0301 }
0302
0303 typedef edm::PSimHitContainer::const_iterator SimHitIter;
0304 typedef RPCRecHitCollection::const_iterator RecHitIter;
0305 typedef std::vector<TrackPSimHitRef> SimHitRefs;
0306
0307
0308 SimHitRefs muonSimHits;
0309
0310 for (int i = 0, n = simParticleHandle->size(); i < n; ++i) {
0311 TrackingParticleRef simParticle(simParticleHandle, i);
0312 if (simParticle->pt() < 1.0 or simParticle->p() < 2.5)
0313 continue;
0314
0315
0316 SimHitRefs simHitsFromParticle;
0317 auto range = std::equal_range(simHitsTPAssoc->begin(),
0318 simHitsTPAssoc->end(),
0319 std::make_pair(simParticle, TrackPSimHitRef()),
0320 SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0321 for (auto simParticleToHit = range.first; simParticleToHit != range.second; ++simParticleToHit) {
0322 auto simHit = simParticleToHit->second;
0323 const DetId detId(simHit->detUnitId());
0324 if (detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC)
0325 continue;
0326
0327 simHitsFromParticle.push_back(simParticleToHit->second);
0328 }
0329 const int nRPCHit = simHitsFromParticle.size();
0330 const bool hasRPCHit = nRPCHit > 0;
0331
0332 if (abs(simParticle->pdgId()) == 13) {
0333 muonSimHits.insert(muonSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
0334
0335
0336 int nRPCHitBarrel = 0;
0337 int nRPCHitEndcap = 0;
0338 for (const auto &simHit : simHitsFromParticle) {
0339 const RPCDetId rpcDetId{simHit->detUnitId()};
0340 const RPCRoll *roll = rpcGeom->roll(rpcDetId);
0341 if (!roll)
0342 continue;
0343
0344 if (rpcDetId.region() == 0)
0345 ++nRPCHitBarrel;
0346 else
0347 ++nRPCHitEndcap;
0348 }
0349
0350
0351 h_nRPCHitPerSimMuon->Fill(nRPCHit);
0352 if (nRPCHitBarrel and nRPCHitEndcap) {
0353 h_nRPCHitPerSimMuonOverlap->Fill(nRPCHit);
0354 h_simMuonOverlap_pt->Fill(simParticle->pt());
0355 h_simMuonOverlap_eta->Fill(simParticle->eta());
0356 h_simMuonOverlap_phi->Fill(simParticle->phi());
0357 } else if (nRPCHitBarrel) {
0358 h_nRPCHitPerSimMuonBarrel->Fill(nRPCHit);
0359 h_simMuonBarrel_pt->Fill(simParticle->pt());
0360 h_simMuonBarrel_eta->Fill(simParticle->eta());
0361 h_simMuonBarrel_phi->Fill(simParticle->phi());
0362 } else if (nRPCHitEndcap) {
0363 h_nRPCHitPerSimMuonEndcap->Fill(nRPCHit);
0364 h_simMuonEndcap_pt->Fill(simParticle->pt());
0365 h_simMuonEndcap_eta->Fill(simParticle->eta());
0366 h_simMuonEndcap_phi->Fill(simParticle->phi());
0367 } else {
0368 h_simMuonNoRPC_pt->Fill(simParticle->pt());
0369 h_simMuonNoRPC_eta->Fill(simParticle->eta());
0370 h_simMuonNoRPC_phi->Fill(simParticle->phi());
0371 }
0372 }
0373
0374 if (hasRPCHit) {
0375 switch (simParticle->pdgId()) {
0376 case 13:
0377 h_simParticleType->Fill(0);
0378 break;
0379 case -13:
0380 h_simParticleType->Fill(1);
0381 break;
0382 case 11:
0383 h_simParticleType->Fill(2);
0384 break;
0385 case -11:
0386 h_simParticleType->Fill(3);
0387 break;
0388 case 211:
0389 h_simParticleType->Fill(4);
0390 break;
0391 case -211:
0392 h_simParticleType->Fill(5);
0393 break;
0394 case 321:
0395 h_simParticleType->Fill(6);
0396 break;
0397 case -321:
0398 h_simParticleType->Fill(7);
0399 break;
0400 case 2212:
0401 h_simParticleType->Fill(8);
0402 break;
0403 case -2212:
0404 h_simParticleType->Fill(9);
0405 break;
0406 default:
0407 h_simParticleType->Fill(10);
0408 break;
0409 }
0410 }
0411 }
0412
0413
0414 int nRefHitBarrel = 0, nRefHitEndcap = 0;
0415 for (const auto &simHit : muonSimHits) {
0416 const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
0417 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
0418
0419 const int region = roll->id().region();
0420 const int ring = roll->id().ring();
0421 const int station = roll->id().station();
0422
0423 if (region == 0) {
0424 ++nRefHitBarrel;
0425 h_.refHitOccupancyBarrel_wheel->Fill(ring);
0426 h_.refHitOccupancyBarrel_station->Fill(station);
0427 h_.refHitOccupancyBarrel_wheel_station->Fill(ring, station);
0428
0429 h_refOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[simHit->detUnitId()]);
0430 } else {
0431 ++nRefHitEndcap;
0432 h_.refHitOccupancyEndcap_disk->Fill(region * station);
0433 h_.refHitOccupancyEndcap_disk_ring->Fill(region * station, ring);
0434
0435 h_refOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[simHit->detUnitId()]);
0436 }
0437 }
0438
0439 h_.nRefHitBarrel->Fill(nRefHitBarrel);
0440 h_.nRefHitEndcap->Fill(nRefHitEndcap);
0441
0442
0443 int sumClusterSizeBarrel = 0, sumClusterSizeEndcap = 0;
0444 int nRecHitBarrel = 0, nRecHitEndcap = 0;
0445 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0446 const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
0447 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
0448 if (!roll)
0449 continue;
0450
0451 const int region = roll->id().region();
0452 const int ring = roll->id().ring();
0453 const int station = roll->id().station();
0454
0455 const double time = recHitIter->timeError() >= 0 ? recHitIter->time() : recHitIter->BunchX() * 25;
0456
0457 h_.clusterSize->Fill(recHitIter->clusterSize());
0458
0459 if (region == 0) {
0460 ++nRecHitBarrel;
0461 sumClusterSizeBarrel += recHitIter->clusterSize();
0462 h_.clusterSizeBarrel->Fill(recHitIter->clusterSize());
0463 h_.recHitOccupancyBarrel_wheel->Fill(ring);
0464 h_.recHitOccupancyBarrel_station->Fill(station);
0465 h_.recHitOccupancyBarrel_wheel_station->Fill(ring, station);
0466
0467 h_allOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[detId.rawId()]);
0468
0469 h_.timeBarrel->Fill(time);
0470 } else {
0471 ++nRecHitEndcap;
0472 sumClusterSizeEndcap += recHitIter->clusterSize();
0473 h_.clusterSizeEndcap->Fill(recHitIter->clusterSize());
0474 h_.recHitOccupancyEndcap_disk->Fill(region * station);
0475 h_.recHitOccupancyEndcap_disk_ring->Fill(region * station, ring);
0476
0477 h_allOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[detId.rawId()]);
0478
0479 h_.timeEndcap->Fill(time);
0480 }
0481
0482 if (roll->isIRPC()) {
0483 h_.timeIRPC->Fill(time);
0484 } else {
0485 h_.timeCRPC->Fill(time);
0486 }
0487 }
0488 const double nRecHit = nRecHitBarrel + nRecHitEndcap;
0489 h_.nRecHitBarrel->Fill(nRecHitBarrel);
0490 h_.nRecHitEndcap->Fill(nRecHitEndcap);
0491 if (nRecHit > 0) {
0492 const int sumClusterSize = sumClusterSizeBarrel + sumClusterSizeEndcap;
0493 h_.avgClusterSize->Fill(double(sumClusterSize) / nRecHit);
0494
0495 if (nRecHitBarrel > 0) {
0496 h_.avgClusterSizeBarrel->Fill(double(sumClusterSizeBarrel) / nRecHitBarrel);
0497 }
0498 if (nRecHitEndcap > 0) {
0499 h_.avgClusterSizeEndcap->Fill(double(sumClusterSizeEndcap) / nRecHitEndcap);
0500 }
0501 }
0502
0503
0504 typedef std::map<TrackPSimHitRef, RecHitIter> SimToRecHitMap;
0505 SimToRecHitMap simToRecHitMap;
0506
0507 for (const auto &simHit : muonSimHits) {
0508 const RPCDetId simDetId{simHit->detUnitId()};
0509
0510 const double simX = simHit->localPosition().x();
0511
0512 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0513 const RPCDetId recDetId{recHitIter->rpcId()};
0514 const RPCRoll *recRoll = rpcGeom->roll(recDetId);
0515 if (!recRoll)
0516 continue;
0517
0518 if (simDetId != recDetId)
0519 continue;
0520
0521 const double recX = recHitIter->localPosition().x();
0522 const double newDx = fabs(recX - simX);
0523
0524
0525 SimToRecHitMap::const_iterator prevSimToReco = simToRecHitMap.find(simHit);
0526 if (prevSimToReco == simToRecHitMap.end()) {
0527 simToRecHitMap.insert(std::make_pair(simHit, recHitIter));
0528 } else {
0529 const double oldDx = fabs(prevSimToReco->second->localPosition().x() - simX);
0530
0531 if (newDx < oldDx) {
0532 simToRecHitMap[simHit] = recHitIter;
0533 }
0534 }
0535 }
0536 }
0537
0538
0539
0540 int nMatchHitBarrel = 0, nMatchHitEndcap = 0;
0541 for (const auto &match : simToRecHitMap) {
0542 TrackPSimHitRef simHit = match.first;
0543 RecHitIter recHitIter = match.second;
0544
0545 const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
0546 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
0547
0548 const int region = roll->id().region();
0549 const int ring = roll->id().ring();
0550 const int station = roll->id().station();
0551
0552 const double simX = simHit->localPosition().x();
0553 const double recX = recHitIter->localPosition().x();
0554 const double errX = sqrt(recHitIter->localPositionError().xx());
0555 const double dX = recX - simX;
0556 const double pull = errX == 0 ? -999 : dX / errX;
0557
0558 if (region == 0) {
0559 ++nMatchHitBarrel;
0560 h_.resBarrel->Fill(dX);
0561 h_.pullBarrel->Fill(pull);
0562 h_.matchOccupancyBarrel_wheel->Fill(ring);
0563 h_.matchOccupancyBarrel_station->Fill(station);
0564 h_.matchOccupancyBarrel_wheel_station->Fill(ring, station);
0565
0566 h_.res_wheel_res->Fill(ring, dX);
0567 h_.res_station_res->Fill(station, dX);
0568 h_.pull_wheel_pull->Fill(ring, pull);
0569 h_.pull_station_pull->Fill(station, pull);
0570
0571 h_matchOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[detId.rawId()]);
0572 } else {
0573 ++nMatchHitEndcap;
0574 h_.resEndcap->Fill(dX);
0575 h_.pullEndcap->Fill(pull);
0576 h_.matchOccupancyEndcap_disk->Fill(region * station);
0577 h_.matchOccupancyEndcap_disk_ring->Fill(region * station, ring);
0578
0579 h_.res_disk_res->Fill(region * station, dX);
0580 h_.res_ring_res->Fill(ring, dX);
0581 h_.pull_disk_pull->Fill(region * station, pull);
0582 h_.pull_ring_pull->Fill(ring, pull);
0583
0584 h_matchOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[detId.rawId()]);
0585 }
0586 }
0587 h_.nMatchHitBarrel->Fill(nMatchHitBarrel);
0588 h_.nMatchHitEndcap->Fill(nMatchHitEndcap);
0589
0590
0591 for (reco::MuonCollection::const_iterator muon = muonHandle->begin(); muon != muonHandle->end(); ++muon) {
0592 if (!muon->isGlobalMuon())
0593 continue;
0594
0595 int nRPCHitBarrel = 0;
0596 int nRPCHitEndcap = 0;
0597
0598 const reco::TrackRef glbTrack = muon->globalTrack();
0599 for (trackingRecHit_iterator recHit = glbTrack->recHitsBegin(); recHit != glbTrack->recHitsEnd(); ++recHit) {
0600 if (!(*recHit)->isValid())
0601 continue;
0602 const DetId detId = (*recHit)->geographicalId();
0603 if (detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC)
0604 continue;
0605 const RPCDetId rpcDetId = static_cast<const RPCDetId>(detId);
0606
0607 if (rpcDetId.region() == 0)
0608 ++nRPCHitBarrel;
0609 else
0610 ++nRPCHitEndcap;
0611 }
0612
0613 const int nRPCHit = nRPCHitBarrel + nRPCHitEndcap;
0614 h_nRPCHitPerRecoMuon->Fill(nRPCHit);
0615 if (nRPCHitBarrel and nRPCHitEndcap) {
0616 h_nRPCHitPerRecoMuonOverlap->Fill(nRPCHit);
0617 h_recoMuonOverlap_pt->Fill(muon->pt());
0618 h_recoMuonOverlap_eta->Fill(muon->eta());
0619 h_recoMuonOverlap_phi->Fill(muon->phi());
0620 } else if (nRPCHitBarrel) {
0621 h_nRPCHitPerRecoMuonBarrel->Fill(nRPCHit);
0622 h_recoMuonBarrel_pt->Fill(muon->pt());
0623 h_recoMuonBarrel_eta->Fill(muon->eta());
0624 h_recoMuonBarrel_phi->Fill(muon->phi());
0625 } else if (nRPCHitEndcap) {
0626 h_nRPCHitPerRecoMuonEndcap->Fill(nRPCHit);
0627 h_recoMuonEndcap_pt->Fill(muon->pt());
0628 h_recoMuonEndcap_eta->Fill(muon->eta());
0629 h_recoMuonEndcap_phi->Fill(muon->phi());
0630 } else {
0631 h_recoMuonNoRPC_pt->Fill(muon->pt());
0632 h_recoMuonNoRPC_eta->Fill(muon->eta());
0633 h_recoMuonNoRPC_phi->Fill(muon->phi());
0634 }
0635 }
0636
0637 h_eventCount->Fill(2);
0638 }
0639
0640 DEFINE_FWK_MODULE(RPCRecHitValid);
