File indexing completed on 2023-10-25 10:07:36
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 h_refPunchOccupancyBarrel_wheel =
0170 booker.book1D("RefPunchOccupancyBarrel_wheel", "RefPunchthrough occupancy", 5, -2.5, 2.5);
0171 h_refPunchOccupancyEndcap_disk =
0172 booker.book1D("RefPunchOccupancyEndcap_disk", "RefPunchthrough occupancy", 9, -4.5, 4.5);
0173 h_refPunchOccupancyBarrel_station =
0174 booker.book1D("RefPunchOccupancyBarrel_station", "RefPunchthrough occupancy", 4, 0.5, 4.5);
0175 h_recPunchOccupancyBarrel_wheel =
0176 booker.book1D("RecPunchOccupancyBarrel_wheel", "Punchthrough recHit occupancy", 5, -2.5, 2.5);
0177 h_recPunchOccupancyEndcap_disk =
0178 booker.book1D("RecPunchOccupancyEndcap_disk", "Punchthrough recHit occupancy", 9, -4.5, 4.5);
0179 h_recPunchOccupancyBarrel_station =
0180 booker.book1D("RecPunchOccupancyBarrel_station", "Punchthrough recHit occupancy", 4, 0.5, 4.5);
0181
0182 h_refPunchOccupancyBarrel_wheel->getTH1()->SetMinimum(0);
0183 h_refPunchOccupancyEndcap_disk->getTH1()->SetMinimum(0);
0184 h_refPunchOccupancyBarrel_station->getTH1()->SetMinimum(0);
0185 h_recPunchOccupancyBarrel_wheel->getTH1()->SetMinimum(0);
0186 h_recPunchOccupancyEndcap_disk->getTH1()->SetMinimum(0);
0187 h_recPunchOccupancyBarrel_station->getTH1()->SetMinimum(0);
0188
0189 h_refPunchOccupancyBarrel_wheel_station =
0190 booker.book2D("RefPunchOccupancyBarrel_wheel_station", "RefPunchthrough occupancy", 5, -2.5, 2.5, 4, 0.5, 4.5);
0191 h_refPunchOccupancyEndcap_disk_ring =
0192 booker.book2D("RefPunchOccupancyEndcap_disk_ring", "RefPunchthrough occupancy", 9, -4.5, 4.5, 4, 0.5, 4.5);
0193 h_recPunchOccupancyBarrel_wheel_station = booker.book2D(
0194 "RecPunchOccupancyBarrel_wheel_station", "Punchthrough recHit occupancy", 5, -2.5, 2.5, 4, 0.5, 4.5);
0195 h_recPunchOccupancyEndcap_disk_ring =
0196 booker.book2D("RecPunchOccupancyEndcap_disk_ring", "Punchthrough recHit occupancy", 9, -4.5, 4.5, 4, 0.5, 4.5);
0197
0198 h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetOption("COLZ");
0199 h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetOption("COLZ");
0200 h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetOption("COLZ");
0201 h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetOption("COLZ");
0202
0203 h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetContour(10);
0204 h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetContour(10);
0205 h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetContour(10);
0206 h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetContour(10);
0207
0208 h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetStats(false);
0209 h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetStats(false);
0210 h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetStats(false);
0211 h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetStats(false);
0212
0213 h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetMinimum(0);
0214 h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetMinimum(0);
0215 h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetMinimum(0);
0216 h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetMinimum(0);
0217
0218 for (int i = 1; i <= 5; ++i) {
0219 TString binLabel = Form("Wheel %d", i - 3);
0220 h_refPunchOccupancyBarrel_wheel->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
0221 h_refPunchOccupancyBarrel_wheel_station->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
0222 h_recPunchOccupancyBarrel_wheel->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
0223 h_recPunchOccupancyBarrel_wheel_station->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
0224 }
0225
0226 for (int i = 1; i <= 9; ++i) {
0227 TString binLabel = Form("Disk %d", i - 5);
0228 h_refPunchOccupancyEndcap_disk->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
0229 h_refPunchOccupancyEndcap_disk_ring->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
0230 h_recPunchOccupancyEndcap_disk->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
0231 h_recPunchOccupancyEndcap_disk_ring->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
0232 }
0233
0234 for (int i = 1; i <= 4; ++i) {
0235 TString binLabel = Form("Station %d", i);
0236 h_refPunchOccupancyBarrel_station->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
0237 h_refPunchOccupancyBarrel_wheel_station->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
0238 h_recPunchOccupancyBarrel_station->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
0239 h_recPunchOccupancyBarrel_wheel_station->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
0240 }
0241
0242 for (int i = 1; i <= 4; ++i) {
0243 TString binLabel = Form("Ring %d", i);
0244 h_refPunchOccupancyEndcap_disk_ring->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
0245 h_recPunchOccupancyEndcap_disk_ring->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
0246 }
0247
0248
0249 auto rpcGeom = eventSetup.getHandle(rpcGeomTokenInRun_);
0250
0251 int nRPCRollBarrel = 0, nRPCRollEndcap = 0;
0252
0253 TrackingGeometry::DetContainer rpcDets = rpcGeom->dets();
0254 for (auto det : rpcDets) {
0255 auto rpcCh = dynamic_cast<const RPCChamber *>(det);
0256 if (!rpcCh)
0257 continue;
0258
0259 std::vector<const RPCRoll *> rolls = rpcCh->rolls();
0260 for (auto roll : rolls) {
0261 if (!roll)
0262 continue;
0263
0264
0265 const int rawId = roll->geographicalId().rawId();
0266
0267
0268
0269 if (roll->isBarrel()) {
0270 detIdToIndexMapBarrel_[rawId] = nRPCRollBarrel;
0271
0272 ++nRPCRollBarrel;
0273 } else {
0274 detIdToIndexMapEndcap_[rawId] = nRPCRollEndcap;
0275
0276 ++nRPCRollEndcap;
0277 }
0278 }
0279 }
0280
0281 booker.setCurrentFolder(subDir_ + "/Occupancy");
0282 h_matchOccupancyBarrel_detId = booker.book1D("MatchOccupancyBarrel_detId",
0283 "Matched hit occupancy;roll index (can be arbitrary)",
0284 nRPCRollBarrel,
0285 0,
0286 nRPCRollBarrel);
0287 h_matchOccupancyEndcap_detId = booker.book1D("MatchOccupancyEndcap_detId",
0288 "Matched hit occupancy;roll index (can be arbitrary)",
0289 nRPCRollEndcap,
0290 0,
0291 nRPCRollEndcap);
0292 h_refOccupancyBarrel_detId = booker.book1D("RefOccupancyBarrel_detId",
0293 "Reference hit occupancy;roll index (can be arbitrary)",
0294 nRPCRollBarrel,
0295 0,
0296 nRPCRollBarrel);
0297 h_refOccupancyEndcap_detId = booker.book1D("RefOccupancyEndcap_detId",
0298 "Reference hit occupancy;roll index (can be arbitrary)",
0299 nRPCRollEndcap,
0300 0,
0301 nRPCRollEndcap);
0302 h_noiseOccupancyBarrel_detId = booker.book1D(
0303 "NoiseOccupancyBarrel_detId", "Noise occupancy;roll index (can be arbitrary)", nRPCRollBarrel, 0, nRPCRollBarrel);
0304 h_noiseOccupancyEndcap_detId = booker.book1D(
0305 "NoiseOccupancyEndcap_detId", "Noise occupancy;roll index (can be arbitrary)", nRPCRollEndcap, 0, nRPCRollEndcap);
0306
0307 h_matchOccupancyBarrel_detId->getTH1()->SetMinimum(0);
0308 h_matchOccupancyEndcap_detId->getTH1()->SetMinimum(0);
0309 h_refOccupancyBarrel_detId->getTH1()->SetMinimum(0);
0310 h_refOccupancyEndcap_detId->getTH1()->SetMinimum(0);
0311 h_noiseOccupancyBarrel_detId->getTH1()->SetMinimum(0);
0312 h_noiseOccupancyEndcap_detId->getTH1()->SetMinimum(0);
0313
0314 h_rollAreaBarrel_detId = booker.bookProfile(
0315 "RollAreaBarrel_detId", "Roll area;roll index;Area", nRPCRollBarrel, 0., 1. * nRPCRollBarrel, 0., 1e5);
0316 h_rollAreaEndcap_detId = booker.bookProfile(
0317 "RollAreaEndcap_detId", "Roll area;roll index;Area", nRPCRollEndcap, 0., 1. * nRPCRollEndcap, 0., 1e5);
0318
0319 for (auto detIdToIndex : detIdToIndexMapBarrel_) {
0320 const int rawId = detIdToIndex.first;
0321 const int index = detIdToIndex.second;
0322
0323 const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
0324 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
0325
0326
0327
0328
0329
0330 const StripTopology &topol = roll->specificTopology();
0331 const double area = topol.stripLength() * topol.nstrips() * topol.pitch();
0332
0333 h_rollAreaBarrel_detId->Fill(index, area);
0334 }
0335
0336 for (auto detIdToIndex : detIdToIndexMapEndcap_) {
0337 const int rawId = detIdToIndex.first;
0338 const int index = detIdToIndex.second;
0339
0340 const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
0341 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
0342
0343
0344
0345
0346
0347 const StripTopology &topol = roll->specificTopology();
0348 const double area = topol.stripLength() * topol.nstrips() * topol.pitch();
0349
0350 h_rollAreaEndcap_detId->Fill(index, area);
0351 }
0352 }
0353
0354 void RPCRecHitValid::analyze(const edm::Event &event, const edm::EventSetup &eventSetup) {
0355 h_eventCount->Fill(1);
0356
0357
0358 auto rpcGeom = eventSetup.getHandle(rpcGeomToken_);
0359
0360
0361 edm::Handle<edm::PSimHitContainer> simHitHandle;
0362 if (!event.getByToken(simHitToken_, simHitHandle)) {
0363 edm::LogInfo("RPCRecHitValid") << "Cannot find simHit collection\n";
0364 return;
0365 }
0366
0367
0368 edm::Handle<RPCRecHitCollection> recHitHandle;
0369 if (!event.getByToken(recHitToken_, recHitHandle)) {
0370 edm::LogInfo("RPCRecHitValid") << "Cannot find recHit collection\n";
0371 return;
0372 }
0373
0374
0375 edm::Handle<TrackingParticleCollection> simParticleHandle;
0376 if (!event.getByToken(simParticleToken_, simParticleHandle)) {
0377 edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle collection\n";
0378 return;
0379 }
0380
0381
0382 edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0383 if (!event.getByToken(simHitAssocToken_, simHitsTPAssoc)) {
0384 edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle to SimHit association map\n";
0385 return;
0386 }
0387
0388
0389 edm::Handle<reco::MuonCollection> muonHandle;
0390 if (!event.getByToken(muonToken_, muonHandle)) {
0391 edm::LogInfo("RPCRecHitValid") << "Cannot find muon collection\n";
0392 return;
0393 }
0394
0395 typedef edm::PSimHitContainer::const_iterator SimHitIter;
0396 typedef RPCRecHitCollection::const_iterator RecHitIter;
0397 typedef std::vector<TrackPSimHitRef> SimHitRefs;
0398
0399
0400 SimHitRefs muonSimHits, pthrSimHits;
0401
0402 for (int i = 0, n = simParticleHandle->size(); i < n; ++i) {
0403 TrackingParticleRef simParticle(simParticleHandle, i);
0404 if (simParticle->pt() < 1.0 or simParticle->p() < 2.5)
0405 continue;
0406
0407
0408 SimHitRefs simHitsFromParticle;
0409 auto range = std::equal_range(simHitsTPAssoc->begin(),
0410 simHitsTPAssoc->end(),
0411 std::make_pair(simParticle, TrackPSimHitRef()),
0412 SimHitTPAssociationProducer::simHitTPAssociationListGreater);
0413 for (auto simParticleToHit = range.first; simParticleToHit != range.second; ++simParticleToHit) {
0414 auto simHit = simParticleToHit->second;
0415 const DetId detId(simHit->detUnitId());
0416 if (detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC)
0417 continue;
0418
0419 simHitsFromParticle.push_back(simParticleToHit->second);
0420 }
0421 const int nRPCHit = simHitsFromParticle.size();
0422 const bool hasRPCHit = nRPCHit > 0;
0423
0424 if (abs(simParticle->pdgId()) == 13) {
0425 muonSimHits.insert(muonSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
0426
0427
0428 int nRPCHitBarrel = 0;
0429 int nRPCHitEndcap = 0;
0430 for (const auto &simHit : simHitsFromParticle) {
0431 const RPCDetId rpcDetId = static_cast<const RPCDetId>(simHit->detUnitId());
0432 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
0433 if (!roll)
0434 continue;
0435
0436 if (rpcDetId.region() == 0)
0437 ++nRPCHitBarrel;
0438 else
0439 ++nRPCHitEndcap;
0440 }
0441
0442
0443 h_nRPCHitPerSimMuon->Fill(nRPCHit);
0444 if (nRPCHitBarrel and nRPCHitEndcap) {
0445 h_nRPCHitPerSimMuonOverlap->Fill(nRPCHit);
0446 h_simMuonOverlap_pt->Fill(simParticle->pt());
0447 h_simMuonOverlap_eta->Fill(simParticle->eta());
0448 h_simMuonOverlap_phi->Fill(simParticle->phi());
0449 } else if (nRPCHitBarrel) {
0450 h_nRPCHitPerSimMuonBarrel->Fill(nRPCHit);
0451 h_simMuonBarrel_pt->Fill(simParticle->pt());
0452 h_simMuonBarrel_eta->Fill(simParticle->eta());
0453 h_simMuonBarrel_phi->Fill(simParticle->phi());
0454 } else if (nRPCHitEndcap) {
0455 h_nRPCHitPerSimMuonEndcap->Fill(nRPCHit);
0456 h_simMuonEndcap_pt->Fill(simParticle->pt());
0457 h_simMuonEndcap_eta->Fill(simParticle->eta());
0458 h_simMuonEndcap_phi->Fill(simParticle->phi());
0459 } else {
0460 h_simMuonNoRPC_pt->Fill(simParticle->pt());
0461 h_simMuonNoRPC_eta->Fill(simParticle->eta());
0462 h_simMuonNoRPC_phi->Fill(simParticle->phi());
0463 }
0464 } else {
0465 pthrSimHits.insert(pthrSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
0466 }
0467
0468 if (hasRPCHit) {
0469 switch (simParticle->pdgId()) {
0470 case 13:
0471 h_simParticleType->Fill(0);
0472 break;
0473 case -13:
0474 h_simParticleType->Fill(1);
0475 break;
0476 case 11:
0477 h_simParticleType->Fill(2);
0478 break;
0479 case -11:
0480 h_simParticleType->Fill(3);
0481 break;
0482 case 211:
0483 h_simParticleType->Fill(4);
0484 break;
0485 case -211:
0486 h_simParticleType->Fill(5);
0487 break;
0488 case 321:
0489 h_simParticleType->Fill(6);
0490 break;
0491 case -321:
0492 h_simParticleType->Fill(7);
0493 break;
0494 case 2212:
0495 h_simParticleType->Fill(8);
0496 break;
0497 case -2212:
0498 h_simParticleType->Fill(9);
0499 break;
0500 default:
0501 h_simParticleType->Fill(10);
0502 break;
0503 }
0504 }
0505 }
0506
0507
0508 int nRefHitBarrel = 0, nRefHitEndcap = 0;
0509 for (const auto &simHit : muonSimHits) {
0510 const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
0511 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
0512
0513 const int region = roll->id().region();
0514 const int ring = roll->id().ring();
0515
0516 const int station = roll->id().station();
0517
0518
0519
0520 if (region == 0) {
0521 ++nRefHitBarrel;
0522 h_.refHitOccupancyBarrel_wheel->Fill(ring);
0523 h_.refHitOccupancyBarrel_station->Fill(station);
0524 h_.refHitOccupancyBarrel_wheel_station->Fill(ring, station);
0525
0526 h_refOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[simHit->detUnitId()]);
0527 } else {
0528 ++nRefHitEndcap;
0529 h_.refHitOccupancyEndcap_disk->Fill(region * station);
0530 h_.refHitOccupancyEndcap_disk_ring->Fill(region * station, ring);
0531
0532 h_refOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[simHit->detUnitId()]);
0533 }
0534 }
0535
0536
0537
0538 for (const auto &simHit : pthrSimHits) {
0539 const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
0540 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
0541
0542 const int region = roll->id().region();
0543 const int ring = roll->id().ring();
0544
0545 const int station = roll->id().station();
0546
0547
0548
0549 if (region == 0) {
0550 ++nRefHitBarrel;
0551 h_refPunchOccupancyBarrel_wheel->Fill(ring);
0552 h_refPunchOccupancyBarrel_station->Fill(station);
0553 h_refPunchOccupancyBarrel_wheel_station->Fill(ring, station);
0554
0555 h_refOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[simHit->detUnitId()]);
0556 } else {
0557 ++nRefHitEndcap;
0558 h_refPunchOccupancyEndcap_disk->Fill(region * station);
0559 h_refPunchOccupancyEndcap_disk_ring->Fill(region * station, ring);
0560
0561 h_refOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[simHit->detUnitId()]);
0562 }
0563 }
0564 h_.nRefHitBarrel->Fill(nRefHitBarrel);
0565 h_.nRefHitEndcap->Fill(nRefHitEndcap);
0566
0567
0568 int sumClusterSizeBarrel = 0, sumClusterSizeEndcap = 0;
0569 int nRecHitBarrel = 0, nRecHitEndcap = 0;
0570 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0571 const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
0572 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
0573 if (!roll)
0574 continue;
0575
0576 const int region = roll->id().region();
0577 const int ring = roll->id().ring();
0578
0579 const int station = roll->id().station();
0580
0581
0582
0583 const double time = recHitIter->timeError() >= 0 ? recHitIter->time() : recHitIter->BunchX() * 25;
0584
0585 h_.clusterSize->Fill(recHitIter->clusterSize());
0586
0587 if (region == 0) {
0588 ++nRecHitBarrel;
0589 sumClusterSizeBarrel += recHitIter->clusterSize();
0590 h_.clusterSizeBarrel->Fill(recHitIter->clusterSize());
0591 h_.recHitOccupancyBarrel_wheel->Fill(ring);
0592 h_.recHitOccupancyBarrel_station->Fill(station);
0593 h_.recHitOccupancyBarrel_wheel_station->Fill(ring, station);
0594
0595 h_.timeBarrel->Fill(time);
0596 } else {
0597 ++nRecHitEndcap;
0598 sumClusterSizeEndcap += recHitIter->clusterSize();
0599 h_.clusterSizeEndcap->Fill(recHitIter->clusterSize());
0600 h_.recHitOccupancyEndcap_disk->Fill(region * station);
0601 h_.recHitOccupancyEndcap_disk_ring->Fill(region * station, ring);
0602
0603 h_.timeEndcap->Fill(time);
0604 }
0605
0606 if (roll->isIRPC()) {
0607 h_.timeIRPC->Fill(time);
0608 } else {
0609 h_.timeCRPC->Fill(time);
0610 }
0611 }
0612 const double nRecHit = nRecHitBarrel + nRecHitEndcap;
0613 h_.nRecHitBarrel->Fill(nRecHitBarrel);
0614 h_.nRecHitEndcap->Fill(nRecHitEndcap);
0615 if (nRecHit > 0) {
0616 const int sumClusterSize = sumClusterSizeBarrel + sumClusterSizeEndcap;
0617 h_.avgClusterSize->Fill(double(sumClusterSize) / nRecHit);
0618
0619 if (nRecHitBarrel > 0) {
0620 h_.avgClusterSizeBarrel->Fill(double(sumClusterSizeBarrel) / nRecHitBarrel);
0621 }
0622 if (nRecHitEndcap > 0) {
0623 h_.avgClusterSizeEndcap->Fill(double(sumClusterSizeEndcap) / nRecHitEndcap);
0624 }
0625 }
0626
0627
0628 typedef std::map<TrackPSimHitRef, RecHitIter> SimToRecHitMap;
0629 SimToRecHitMap simToRecHitMap;
0630
0631 for (const auto &simHit : muonSimHits) {
0632 const RPCDetId simDetId = static_cast<const RPCDetId>(simHit->detUnitId());
0633
0634
0635
0636 const double simX = simHit->localPosition().x();
0637
0638 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0639 const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
0640 const RPCRoll *recRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(recDetId));
0641 if (!recRoll)
0642 continue;
0643
0644 if (simDetId != recDetId)
0645 continue;
0646
0647 const double recX = recHitIter->localPosition().x();
0648 const double newDx = fabs(recX - simX);
0649
0650
0651 SimToRecHitMap::const_iterator prevSimToReco = simToRecHitMap.find(simHit);
0652 if (prevSimToReco == simToRecHitMap.end()) {
0653 simToRecHitMap.insert(std::make_pair(simHit, recHitIter));
0654 } else {
0655 const double oldDx = fabs(prevSimToReco->second->localPosition().x() - simX);
0656
0657 if (newDx < oldDx) {
0658 simToRecHitMap[simHit] = recHitIter;
0659 }
0660 }
0661 }
0662 }
0663
0664
0665
0666 int nMatchHitBarrel = 0, nMatchHitEndcap = 0;
0667 for (const auto &match : simToRecHitMap) {
0668 TrackPSimHitRef simHit = match.first;
0669 RecHitIter recHitIter = match.second;
0670
0671 const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
0672 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
0673
0674 const int region = roll->id().region();
0675 const int ring = roll->id().ring();
0676
0677 const int station = roll->id().station();
0678
0679
0680
0681 const double simX = simHit->localPosition().x();
0682 const double recX = recHitIter->localPosition().x();
0683 const double errX = sqrt(recHitIter->localPositionError().xx());
0684 const double dX = recX - simX;
0685 const double pull = errX == 0 ? -999 : dX / errX;
0686
0687
0688
0689
0690 if (region == 0) {
0691 ++nMatchHitBarrel;
0692 h_.resBarrel->Fill(dX);
0693 h_.pullBarrel->Fill(pull);
0694 h_.matchOccupancyBarrel_wheel->Fill(ring);
0695 h_.matchOccupancyBarrel_station->Fill(station);
0696 h_.matchOccupancyBarrel_wheel_station->Fill(ring, station);
0697
0698 h_.res_wheel_res->Fill(ring, dX);
0699 h_.res_station_res->Fill(station, dX);
0700 h_.pull_wheel_pull->Fill(ring, pull);
0701 h_.pull_station_pull->Fill(station, pull);
0702
0703 h_matchOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[detId.rawId()]);
0704 } else {
0705 ++nMatchHitEndcap;
0706 h_.resEndcap->Fill(dX);
0707 h_.pullEndcap->Fill(pull);
0708 h_.matchOccupancyEndcap_disk->Fill(region * station);
0709 h_.matchOccupancyEndcap_disk_ring->Fill(region * station, ring);
0710
0711 h_.res_disk_res->Fill(region * station, dX);
0712 h_.res_ring_res->Fill(ring, dX);
0713 h_.pull_disk_pull->Fill(region * station, pull);
0714 h_.pull_ring_pull->Fill(ring, pull);
0715
0716 h_matchOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[detId.rawId()]);
0717 }
0718 }
0719 h_.nMatchHitBarrel->Fill(nMatchHitBarrel);
0720 h_.nMatchHitEndcap->Fill(nMatchHitEndcap);
0721
0722
0723 for (reco::MuonCollection::const_iterator muon = muonHandle->begin(); muon != muonHandle->end(); ++muon) {
0724 if (!muon->isGlobalMuon())
0725 continue;
0726
0727 int nRPCHitBarrel = 0;
0728 int nRPCHitEndcap = 0;
0729
0730 const reco::TrackRef glbTrack = muon->globalTrack();
0731 for (trackingRecHit_iterator recHit = glbTrack->recHitsBegin(); recHit != glbTrack->recHitsEnd(); ++recHit) {
0732 if (!(*recHit)->isValid())
0733 continue;
0734 const DetId detId = (*recHit)->geographicalId();
0735 if (detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC)
0736 continue;
0737 const RPCDetId rpcDetId = static_cast<const RPCDetId>(detId);
0738
0739 if (rpcDetId.region() == 0)
0740 ++nRPCHitBarrel;
0741 else
0742 ++nRPCHitEndcap;
0743 }
0744
0745 const int nRPCHit = nRPCHitBarrel + nRPCHitEndcap;
0746 h_nRPCHitPerRecoMuon->Fill(nRPCHit);
0747 if (nRPCHitBarrel and nRPCHitEndcap) {
0748 h_nRPCHitPerRecoMuonOverlap->Fill(nRPCHit);
0749 h_recoMuonOverlap_pt->Fill(muon->pt());
0750 h_recoMuonOverlap_eta->Fill(muon->eta());
0751 h_recoMuonOverlap_phi->Fill(muon->phi());
0752 } else if (nRPCHitBarrel) {
0753 h_nRPCHitPerRecoMuonBarrel->Fill(nRPCHit);
0754 h_recoMuonBarrel_pt->Fill(muon->pt());
0755 h_recoMuonBarrel_eta->Fill(muon->eta());
0756 h_recoMuonBarrel_phi->Fill(muon->phi());
0757 } else if (nRPCHitEndcap) {
0758 h_nRPCHitPerRecoMuonEndcap->Fill(nRPCHit);
0759 h_recoMuonEndcap_pt->Fill(muon->pt());
0760 h_recoMuonEndcap_eta->Fill(muon->eta());
0761 h_recoMuonEndcap_phi->Fill(muon->phi());
0762 } else {
0763 h_recoMuonNoRPC_pt->Fill(muon->pt());
0764 h_recoMuonNoRPC_eta->Fill(muon->eta());
0765 h_recoMuonNoRPC_phi->Fill(muon->phi());
0766 }
0767 }
0768
0769
0770 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0771 const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
0772 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
0773
0774 const int region = roll->id().region();
0775 const int ring = roll->id().ring();
0776
0777 const int station = roll->id().station();
0778
0779
0780
0781 bool matched = false;
0782 for (const auto &match : simToRecHitMap) {
0783 if (recHitIter == match.second) {
0784 matched = true;
0785 break;
0786 }
0787 }
0788
0789 if (!matched) {
0790 int nPunchMatched = 0;
0791
0792 for (const auto &simHit : pthrSimHits) {
0793 const int absSimHitPType = abs(simHit->particleType());
0794 if (absSimHitPType == 13)
0795 continue;
0796
0797 const RPCDetId simDetId = static_cast<const RPCDetId>(simHit->detUnitId());
0798 if (simDetId == detId)
0799 ++nPunchMatched;
0800 }
0801
0802 if (nPunchMatched > 0) {
0803 if (region == 0) {
0804 h_recPunchOccupancyBarrel_wheel->Fill(ring);
0805 h_recPunchOccupancyBarrel_station->Fill(station);
0806 h_recPunchOccupancyBarrel_wheel_station->Fill(ring, station);
0807 } else {
0808 h_recPunchOccupancyEndcap_disk->Fill(region * station);
0809 h_recPunchOccupancyEndcap_disk_ring->Fill(region * station, ring);
0810 }
0811 }
0812 }
0813 }
0814
0815
0816 for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0817 const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
0818 const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(recDetId));
0819
0820 const int region = roll->id().region();
0821
0822
0823
0824
0825
0826
0827 const double recX = recHitIter->localPosition().x();
0828 const double recErrX = sqrt(recHitIter->localPositionError().xx());
0829
0830 bool matched = false;
0831 for (SimHitIter simHitIter = simHitHandle->begin(); simHitIter != simHitHandle->end(); ++simHitIter) {
0832 const RPCDetId simDetId = static_cast<const RPCDetId>(simHitIter->detUnitId());
0833 const RPCRoll *simRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(simDetId));
0834 if (!simRoll)
0835 continue;
0836
0837 if (simDetId != recDetId)
0838 continue;
0839
0840 const double simX = simHitIter->localPosition().x();
0841 const double dX = fabs(recX - simX);
0842
0843 if (dX / recErrX < 5) {
0844 matched = true;
0845 break;
0846 }
0847 }
0848
0849 if (!matched) {
0850 if (region == 0) {
0851 h_noiseOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[recDetId.rawId()]);
0852 } else {
0853 h_noiseOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[recDetId.rawId()]);
0854 }
0855 }
0856 }
0857
0858 h_eventCount->Fill(2);
0859 }
0860
0861 DEFINE_FWK_MODULE(RPCRecHitValid);
