Back to home page

Project CMSSW displayed by LXR

 
 

    


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   // Book MonitorElements
0039   h_.bookHistograms(booker, subDir_);
0040 
0041   // SimHit plots, not compatible to RPCPoint-RPCRecHit comparison
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   // Book roll-by-roll histograms
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       // RPCGeomServ rpcSrv(roll->id());
0265       const int rawId = roll->geographicalId().rawId();
0266       // if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name()
0267       // << ' ' << rawId << endl; continue; }
0268 
0269       if (roll->isBarrel()) {
0270         detIdToIndexMapBarrel_[rawId] = nRPCRollBarrel;
0271         // rollIdToNameMapBarrel_[rawId] = rpcSrv.name();
0272         ++nRPCRollBarrel;
0273       } else {
0274         detIdToIndexMapEndcap_[rawId] = nRPCRollEndcap;
0275         // rollIdToNameMapEndcap_[rawId] = rpcSrv.name();
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     // RPCGeomServ rpcSrv(roll->id());
0327     // if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name() <<
0328     // ' ' << rawId << endl; continue; }
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     // RPCGeomServ rpcSrv(roll->id());
0344     // if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name() <<
0345     // ' ' << rawId << endl; continue; }
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   // Get the RPC Geometry
0358   auto rpcGeom = eventSetup.getHandle(rpcGeomToken_);
0359 
0360   // Retrieve SimHits from the event
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   // Retrieve RecHits from the event
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   // Get SimParticles
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   // Get SimParticle to SimHit association map
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   // Get RecoMuons
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   // TrackingParticles with (and without) RPC simHits
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;  // globalMuon acceptance
0406 
0407     // Collect SimHits from this Tracking Particle
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       // Count number of Barrel hits and Endcap hits
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       // Fill TrackingParticle related histograms
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   // Loop over muon simHits, fill histograms which does not need associations
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     // const int sector = roll->id().sector();
0516     const int station = roll->id().station();
0517     // const int layer = roll->id().layer();
0518     // const int subSector = roll->id().subsector();
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   // Loop over punch-through simHits, fill histograms which does not need
0537   // associations
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     // const int sector = roll->id().sector();
0545     const int station = roll->id().station();
0546     // const int layer = roll->id().layer();
0547     // const int subSector = roll->id().subsector();
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   // Loop over recHits, fill histograms which does not need associations
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     // const int sector = roll->id().sector();
0579     const int station = roll->id().station();
0580     // const int layer = roll->id().layer();
0581     // const int subSector = roll->id().subsector();
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   // Start matching SimHits to RecHits
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     // const RPCRoll* simRoll = dynamic_cast<const
0634     // RPCRoll*>(rpcGeom->roll(simDetId));
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       // Associate SimHit to RecHit
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   // Now we have simHit-recHit mapping
0665   // So we can fill up relavant histograms
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     // const int sector = roll->id().sector();
0677     const int station = roll->id().station();
0678     // const int layer = roll->id().layer();
0679     // const int subsector = roll->id().subsector();
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     // const GlobalPoint simPos = roll->toGlobal(simHitIter->localPosition());
0688     // const GlobalPoint recPos = roll->toGlobal(recHitIter->localPosition());
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   // Reco Muon hits
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   // Find Non-muon hits
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     // const int sector = roll->id().sector();
0777     const int station = roll->id().station();
0778     // const int layer = roll->id().layer();
0779     // const int subsector = roll->id().subsector();
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       // Check if this recHit came from non-muon simHit
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   // Find noise recHits : RecHits without SimHit match
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     //    const int ring = roll->id().ring(); // UNUSED VARIABLE
0822     // const int sector = roll->id().sector();
0823     //    const int station = roll->id().station(); // UNUSED VARIABLE
0824     // const int layer = roll->id().layer();
0825     // const int subsector = roll->id().subsector();
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);