CMS 3D CMS Logo

RPCRecHitValid.cc
Go to the documentation of this file.
3 
6 
17 
18 #include <algorithm>
19 
20 using namespace std;
21 
23 
25  simHitToken_ = consumes<SimHits>(pset.getParameter<edm::InputTag>("simHit"));
26  recHitToken_ = consumes<RecHits>(pset.getParameter<edm::InputTag>("recHit"));
27  simParticleToken_ = consumes<SimParticles>(pset.getParameter<edm::InputTag>("simTrack"));
28  simHitAssocToken_ = consumes<SimHitAssoc>(pset.getParameter<edm::InputTag>("simHitAssoc"));
29  muonToken_ = consumes<reco::MuonCollection>(pset.getParameter<edm::InputTag>("muon"));
30 
31  subDir_ = pset.getParameter<std::string>("subDir");
32 
33  rpcGeomToken_ = esConsumes();
34  rpcGeomTokenInRun_ = esConsumes<edm::Transition::BeginRun>();
35 }
36 
38  // Book MonitorElements
39  h_.bookHistograms(booker, subDir_);
40 
41  // SimHit plots, not compatible to RPCPoint-RPCRecHit comparison
42  booker.setCurrentFolder(subDir_ + "/HitProperty");
43  h_simParticleType = booker.book1D("SimHitPType", "SimHit particle type", 11, 0, 11);
44  h_simParticleType->getTH1()->SetMinimum(0);
45  if (TH1 *h = h_simParticleType->getTH1()) {
46  h->GetXaxis()->SetBinLabel(1, "#mu^{-}");
47  h->GetXaxis()->SetBinLabel(2, "#mu^{+}");
48  h->GetXaxis()->SetBinLabel(3, "e^{-}");
49  h->GetXaxis()->SetBinLabel(4, "e^{+}");
50  h->GetXaxis()->SetBinLabel(5, "#pi^{+}");
51  h->GetXaxis()->SetBinLabel(6, "#pi^{-}");
52  h->GetXaxis()->SetBinLabel(7, "K^{+}");
53  h->GetXaxis()->SetBinLabel(8, "K^{-}");
54  h->GetXaxis()->SetBinLabel(9, "p^{+}");
55  h->GetXaxis()->SetBinLabel(10, "p^{-}");
56  h->GetXaxis()->SetBinLabel(11, "Other");
57  }
58 
59  booker.setCurrentFolder(subDir_ + "/Track");
60 
61  h_nRPCHitPerSimMuon = booker.book1D("NRPCHitPerSimMuon", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
62  h_nRPCHitPerSimMuonBarrel =
63  booker.book1D("NRPCHitPerSimMuonBarrel", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
64  h_nRPCHitPerSimMuonOverlap =
65  booker.book1D("NRPCHitPerSimMuonOverlap", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
66  h_nRPCHitPerSimMuonEndcap =
67  booker.book1D("NRPCHitPerSimMuonEndcap", "Number of RPC SimHit per SimMuon", 11, -0.5, 10.5);
68 
69  h_nRPCHitPerRecoMuon = booker.book1D("NRPCHitPerRecoMuon", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
70  h_nRPCHitPerRecoMuonBarrel =
71  booker.book1D("NRPCHitPerRecoMuonBarrel", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
72  h_nRPCHitPerRecoMuonOverlap =
73  booker.book1D("NRPCHitPerRecoMuonOverlap", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
74  h_nRPCHitPerRecoMuonEndcap =
75  booker.book1D("NRPCHitPerRecoMuonEndcap", "Number of RPC RecHit per RecoMuon", 11, -0.5, 10.5);
76 
77  h_nRPCHitPerSimMuon->getTH1()->SetMinimum(0);
78  h_nRPCHitPerSimMuonBarrel->getTH1()->SetMinimum(0);
79  h_nRPCHitPerSimMuonOverlap->getTH1()->SetMinimum(0);
80  h_nRPCHitPerSimMuonEndcap->getTH1()->SetMinimum(0);
81 
82  h_nRPCHitPerRecoMuon->getTH1()->SetMinimum(0);
83  h_nRPCHitPerRecoMuonBarrel->getTH1()->SetMinimum(0);
84  h_nRPCHitPerRecoMuonOverlap->getTH1()->SetMinimum(0);
85  h_nRPCHitPerRecoMuonEndcap->getTH1()->SetMinimum(0);
86 
87  float ptBins[] = {0, 1, 2, 5, 10, 20, 30, 50, 100, 200, 300, 500};
88  const int nPtBins = sizeof(ptBins) / sizeof(float) - 1;
89  h_simMuonBarrel_pt =
90  booker.book1D("SimMuonBarrel_pt", "SimMuon RPCHit in Barrel p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
91  h_simMuonOverlap_pt =
92  booker.book1D("SimMuonOverlap_pt", "SimMuon RPCHit in Overlap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
93  h_simMuonEndcap_pt =
94  booker.book1D("SimMuonEndcap_pt", "SimMuon RPCHit in Endcap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
95  h_simMuonNoRPC_pt =
96  booker.book1D("SimMuonNoRPC_pt", "SimMuon without RPCHit p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
97  h_simMuonBarrel_eta = booker.book1D("SimMuonBarrel_eta", "SimMuon RPCHit in Barrel #eta;#eta", 50, -2.5, 2.5);
98  h_simMuonOverlap_eta = booker.book1D("SimMuonOverlap_eta", "SimMuon RPCHit in Overlap #eta;#eta", 50, -2.5, 2.5);
99  h_simMuonEndcap_eta = booker.book1D("SimMuonEndcap_eta", "SimMuon RPCHit in Endcap #eta;#eta", 50, -2.5, 2.5);
100  h_simMuonNoRPC_eta = booker.book1D("SimMuonNoRPC_eta", "SimMuon without RPCHit #eta;#eta", 50, -2.5, 2.5);
101  h_simMuonBarrel_phi =
102  booker.book1D("SimMuonBarrel_phi", "SimMuon RPCHit in Barrel #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
103  h_simMuonOverlap_phi =
104  booker.book1D("SimMuonOverlap_phi", "SimMuon RPCHit in Overlap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
105  h_simMuonEndcap_phi =
106  booker.book1D("SimMuonEndcap_phi", "SimMuon RPCHit in Endcap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
107  h_simMuonNoRPC_phi =
108  booker.book1D("SimMuonNoRPC_phi", "SimMuon without RPCHit #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
109 
110  h_recoMuonBarrel_pt =
111  booker.book1D("RecoMuonBarrel_pt", "RecoMuon RPCHit in Barrel p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
112  h_recoMuonOverlap_pt =
113  booker.book1D("RecoMuonOverlap_pt", "RecoMuon RPCHit in Overlap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
114  h_recoMuonEndcap_pt =
115  booker.book1D("RecoMuonEndcap_pt", "RecoMuon RPCHit in Endcap p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
116  h_recoMuonNoRPC_pt =
117  booker.book1D("RecoMuonNoRPC_pt", "RecoMuon without RPCHit p_{T};p_{T} [GeV/c^{2}]", nPtBins, ptBins);
118  h_recoMuonBarrel_eta = booker.book1D("RecoMuonBarrel_eta", "RecoMuon RPCHit in Barrel #eta;#eta", 50, -2.5, 2.5);
119  h_recoMuonOverlap_eta = booker.book1D("RecoMuonOverlap_eta", "RecoMuon RPCHit in Overlap #eta;#eta", 50, -2.5, 2.5);
120  h_recoMuonEndcap_eta = booker.book1D("RecoMuonEndcap_eta", "RecoMuon RPCHit in Endcap #eta;#eta", 50, -2.5, 2.5);
121  h_recoMuonNoRPC_eta = booker.book1D("RecoMuonNoRPC_eta", "RecoMuon without RPCHit #eta;#eta", 50, -2.5, 2.5);
122  h_recoMuonBarrel_phi =
123  booker.book1D("RecoMuonBarrel_phi", "RecoMuon RPCHit in Barrel #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
124  h_recoMuonOverlap_phi =
125  booker.book1D("RecoMuonOverlap_phi", "RecoMuon RPCHit in Overlap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
126  h_recoMuonEndcap_phi =
127  booker.book1D("RecoMuonEndcap_phi", "RecoMuon RPCHit in Endcap #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
128  h_recoMuonNoRPC_phi =
129  booker.book1D("RecoMuonNoRPC_phi", "RecoMuon without RPCHit #phi;#phi", 36, -TMath::Pi(), TMath::Pi());
130 
131  h_simMuonBarrel_pt->getTH1()->SetMinimum(0);
132  h_simMuonOverlap_pt->getTH1()->SetMinimum(0);
133  h_simMuonEndcap_pt->getTH1()->SetMinimum(0);
134  h_simMuonNoRPC_pt->getTH1()->SetMinimum(0);
135  h_simMuonBarrel_eta->getTH1()->SetMinimum(0);
136  h_simMuonOverlap_eta->getTH1()->SetMinimum(0);
137  h_simMuonEndcap_eta->getTH1()->SetMinimum(0);
138  h_simMuonNoRPC_eta->getTH1()->SetMinimum(0);
139  h_simMuonBarrel_phi->getTH1()->SetMinimum(0);
140  h_simMuonOverlap_phi->getTH1()->SetMinimum(0);
141  h_simMuonEndcap_phi->getTH1()->SetMinimum(0);
142  h_simMuonNoRPC_phi->getTH1()->SetMinimum(0);
143 
144  h_recoMuonBarrel_pt->getTH1()->SetMinimum(0);
145  h_recoMuonOverlap_pt->getTH1()->SetMinimum(0);
146  h_recoMuonEndcap_pt->getTH1()->SetMinimum(0);
147  h_recoMuonNoRPC_pt->getTH1()->SetMinimum(0);
148  h_recoMuonBarrel_eta->getTH1()->SetMinimum(0);
149  h_recoMuonOverlap_eta->getTH1()->SetMinimum(0);
150  h_recoMuonEndcap_eta->getTH1()->SetMinimum(0);
151  h_recoMuonNoRPC_eta->getTH1()->SetMinimum(0);
152  h_recoMuonBarrel_phi->getTH1()->SetMinimum(0);
153  h_recoMuonOverlap_phi->getTH1()->SetMinimum(0);
154  h_recoMuonEndcap_phi->getTH1()->SetMinimum(0);
155  h_recoMuonNoRPC_phi->getTH1()->SetMinimum(0);
156 
157  booker.setCurrentFolder(subDir_ + "/Occupancy");
158 
159  h_eventCount = booker.book1D("EventCount", "Event count", 3, 1, 4);
160  h_eventCount->getTH1()->SetMinimum(0);
161  if (h_eventCount) {
162  TH1 *h = h_eventCount->getTH1();
163  h->GetXaxis()->SetBinLabel(1, "eventBegin");
164  h->GetXaxis()->SetBinLabel(2, "eventEnd");
165  h->GetXaxis()->SetBinLabel(3, "run");
166  }
167  h_eventCount->Fill(3);
168 
169  h_refPunchOccupancyBarrel_wheel =
170  booker.book1D("RefPunchOccupancyBarrel_wheel", "RefPunchthrough occupancy", 5, -2.5, 2.5);
171  h_refPunchOccupancyEndcap_disk =
172  booker.book1D("RefPunchOccupancyEndcap_disk", "RefPunchthrough occupancy", 9, -4.5, 4.5);
173  h_refPunchOccupancyBarrel_station =
174  booker.book1D("RefPunchOccupancyBarrel_station", "RefPunchthrough occupancy", 4, 0.5, 4.5);
175  h_recPunchOccupancyBarrel_wheel =
176  booker.book1D("RecPunchOccupancyBarrel_wheel", "Punchthrough recHit occupancy", 5, -2.5, 2.5);
177  h_recPunchOccupancyEndcap_disk =
178  booker.book1D("RecPunchOccupancyEndcap_disk", "Punchthrough recHit occupancy", 9, -4.5, 4.5);
179  h_recPunchOccupancyBarrel_station =
180  booker.book1D("RecPunchOccupancyBarrel_station", "Punchthrough recHit occupancy", 4, 0.5, 4.5);
181 
182  h_refPunchOccupancyBarrel_wheel->getTH1()->SetMinimum(0);
183  h_refPunchOccupancyEndcap_disk->getTH1()->SetMinimum(0);
184  h_refPunchOccupancyBarrel_station->getTH1()->SetMinimum(0);
185  h_recPunchOccupancyBarrel_wheel->getTH1()->SetMinimum(0);
186  h_recPunchOccupancyEndcap_disk->getTH1()->SetMinimum(0);
187  h_recPunchOccupancyBarrel_station->getTH1()->SetMinimum(0);
188 
189  h_refPunchOccupancyBarrel_wheel_station =
190  booker.book2D("RefPunchOccupancyBarrel_wheel_station", "RefPunchthrough occupancy", 5, -2.5, 2.5, 4, 0.5, 4.5);
191  h_refPunchOccupancyEndcap_disk_ring =
192  booker.book2D("RefPunchOccupancyEndcap_disk_ring", "RefPunchthrough occupancy", 9, -4.5, 4.5, 4, 0.5, 4.5);
193  h_recPunchOccupancyBarrel_wheel_station = booker.book2D(
194  "RecPunchOccupancyBarrel_wheel_station", "Punchthrough recHit occupancy", 5, -2.5, 2.5, 4, 0.5, 4.5);
195  h_recPunchOccupancyEndcap_disk_ring =
196  booker.book2D("RecPunchOccupancyEndcap_disk_ring", "Punchthrough recHit occupancy", 9, -4.5, 4.5, 4, 0.5, 4.5);
197 
198  h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetOption("COLZ");
199  h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetOption("COLZ");
200  h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetOption("COLZ");
201  h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetOption("COLZ");
202 
203  h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetContour(10);
204  h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetContour(10);
205  h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetContour(10);
206  h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetContour(10);
207 
208  h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetStats(false);
209  h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetStats(false);
210  h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetStats(false);
211  h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetStats(false);
212 
213  h_refPunchOccupancyBarrel_wheel_station->getTH2F()->SetMinimum(0);
214  h_refPunchOccupancyEndcap_disk_ring->getTH2F()->SetMinimum(0);
215  h_recPunchOccupancyBarrel_wheel_station->getTH2F()->SetMinimum(0);
216  h_recPunchOccupancyEndcap_disk_ring->getTH2F()->SetMinimum(0);
217 
218  for (int i = 1; i <= 5; ++i) {
219  TString binLabel = Form("Wheel %d", i - 3);
220  h_refPunchOccupancyBarrel_wheel->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
221  h_refPunchOccupancyBarrel_wheel_station->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
222  h_recPunchOccupancyBarrel_wheel->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
223  h_recPunchOccupancyBarrel_wheel_station->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
224  }
225 
226  for (int i = 1; i <= 9; ++i) {
227  TString binLabel = Form("Disk %d", i - 5);
228  h_refPunchOccupancyEndcap_disk->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
229  h_refPunchOccupancyEndcap_disk_ring->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
230  h_recPunchOccupancyEndcap_disk->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
231  h_recPunchOccupancyEndcap_disk_ring->getTH2F()->GetXaxis()->SetBinLabel(i, binLabel);
232  }
233 
234  for (int i = 1; i <= 4; ++i) {
235  TString binLabel = Form("Station %d", i);
236  h_refPunchOccupancyBarrel_station->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
237  h_refPunchOccupancyBarrel_wheel_station->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
238  h_recPunchOccupancyBarrel_station->getTH1()->GetXaxis()->SetBinLabel(i, binLabel);
239  h_recPunchOccupancyBarrel_wheel_station->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
240  }
241 
242  for (int i = 1; i <= 4; ++i) {
243  TString binLabel = Form("Ring %d", i);
244  h_refPunchOccupancyEndcap_disk_ring->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
245  h_recPunchOccupancyEndcap_disk_ring->getTH2F()->GetYaxis()->SetBinLabel(i, binLabel);
246  }
247 
248  // Book roll-by-roll histograms
249  auto rpcGeom = eventSetup.getHandle(rpcGeomTokenInRun_);
250 
251  int nRPCRollBarrel = 0, nRPCRollEndcap = 0;
252 
253  TrackingGeometry::DetContainer rpcDets = rpcGeom->dets();
254  for (auto det : rpcDets) {
255  auto rpcCh = dynamic_cast<const RPCChamber *>(det);
256  if (!rpcCh)
257  continue;
258 
259  std::vector<const RPCRoll *> rolls = rpcCh->rolls();
260  for (auto roll : rolls) {
261  if (!roll)
262  continue;
263 
264  // RPCGeomServ rpcSrv(roll->id());
265  const int rawId = roll->geographicalId().rawId();
266  // if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name()
267  // << ' ' << rawId << endl; continue; }
268 
269  if (roll->isBarrel()) {
270  detIdToIndexMapBarrel_[rawId] = nRPCRollBarrel;
271  // rollIdToNameMapBarrel_[rawId] = rpcSrv.name();
272  ++nRPCRollBarrel;
273  } else {
274  detIdToIndexMapEndcap_[rawId] = nRPCRollEndcap;
275  // rollIdToNameMapEndcap_[rawId] = rpcSrv.name();
276  ++nRPCRollEndcap;
277  }
278  }
279  }
280 
281  booker.setCurrentFolder(subDir_ + "/Occupancy");
282  h_matchOccupancyBarrel_detId = booker.book1D("MatchOccupancyBarrel_detId",
283  "Matched hit occupancy;roll index (can be arbitrary)",
284  nRPCRollBarrel,
285  0,
286  nRPCRollBarrel);
287  h_matchOccupancyEndcap_detId = booker.book1D("MatchOccupancyEndcap_detId",
288  "Matched hit occupancy;roll index (can be arbitrary)",
289  nRPCRollEndcap,
290  0,
291  nRPCRollEndcap);
292  h_refOccupancyBarrel_detId = booker.book1D("RefOccupancyBarrel_detId",
293  "Reference hit occupancy;roll index (can be arbitrary)",
294  nRPCRollBarrel,
295  0,
296  nRPCRollBarrel);
297  h_refOccupancyEndcap_detId = booker.book1D("RefOccupancyEndcap_detId",
298  "Reference hit occupancy;roll index (can be arbitrary)",
299  nRPCRollEndcap,
300  0,
301  nRPCRollEndcap);
302  h_noiseOccupancyBarrel_detId = booker.book1D(
303  "NoiseOccupancyBarrel_detId", "Noise occupancy;roll index (can be arbitrary)", nRPCRollBarrel, 0, nRPCRollBarrel);
304  h_noiseOccupancyEndcap_detId = booker.book1D(
305  "NoiseOccupancyEndcap_detId", "Noise occupancy;roll index (can be arbitrary)", nRPCRollEndcap, 0, nRPCRollEndcap);
306 
307  h_matchOccupancyBarrel_detId->getTH1()->SetMinimum(0);
308  h_matchOccupancyEndcap_detId->getTH1()->SetMinimum(0);
309  h_refOccupancyBarrel_detId->getTH1()->SetMinimum(0);
310  h_refOccupancyEndcap_detId->getTH1()->SetMinimum(0);
311  h_noiseOccupancyBarrel_detId->getTH1()->SetMinimum(0);
312  h_noiseOccupancyEndcap_detId->getTH1()->SetMinimum(0);
313 
314  h_rollAreaBarrel_detId = booker.bookProfile(
315  "RollAreaBarrel_detId", "Roll area;roll index;Area", nRPCRollBarrel, 0., 1. * nRPCRollBarrel, 0., 1e5);
316  h_rollAreaEndcap_detId = booker.bookProfile(
317  "RollAreaEndcap_detId", "Roll area;roll index;Area", nRPCRollEndcap, 0., 1. * nRPCRollEndcap, 0., 1e5);
318 
319  for (auto detIdToIndex : detIdToIndexMapBarrel_) {
320  const int rawId = detIdToIndex.first;
321  const int index = detIdToIndex.second;
322 
323  const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
324  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
325 
326  // RPCGeomServ rpcSrv(roll->id());
327  // if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name() <<
328  // ' ' << rawId << endl; continue; }
329 
330  const StripTopology &topol = roll->specificTopology();
331  const double area = topol.stripLength() * topol.nstrips() * topol.pitch();
332 
333  h_rollAreaBarrel_detId->Fill(index, area);
334  }
335 
336  for (auto detIdToIndex : detIdToIndexMapEndcap_) {
337  const int rawId = detIdToIndex.first;
338  const int index = detIdToIndex.second;
339 
340  const RPCDetId rpcDetId = static_cast<const RPCDetId>(rawId);
341  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
342 
343  // RPCGeomServ rpcSrv(roll->id());
344  // if ( !roll->specs()->isRPC() ) { cout << "\nNoRPC : " << rpcSrv.name() <<
345  // ' ' << rawId << endl; continue; }
346 
347  const StripTopology &topol = roll->specificTopology();
348  const double area = topol.stripLength() * topol.nstrips() * topol.pitch();
349 
350  h_rollAreaEndcap_detId->Fill(index, area);
351  }
352 }
353 
355  h_eventCount->Fill(1);
356 
357  // Get the RPC Geometry
358  auto rpcGeom = eventSetup.getHandle(rpcGeomToken_);
359 
360  // Retrieve SimHits from the event
362  if (!event.getByToken(simHitToken_, simHitHandle)) {
363  edm::LogInfo("RPCRecHitValid") << "Cannot find simHit collection\n";
364  return;
365  }
366 
367  // Retrieve RecHits from the event
369  if (!event.getByToken(recHitToken_, recHitHandle)) {
370  edm::LogInfo("RPCRecHitValid") << "Cannot find recHit collection\n";
371  return;
372  }
373 
374  // Get SimParticles
375  edm::Handle<TrackingParticleCollection> simParticleHandle;
376  if (!event.getByToken(simParticleToken_, simParticleHandle)) {
377  edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle collection\n";
378  return;
379  }
380 
381  // Get SimParticle to SimHit association map
383  if (!event.getByToken(simHitAssocToken_, simHitsTPAssoc)) {
384  edm::LogInfo("RPCRecHitValid") << "Cannot find TrackingParticle to SimHit association map\n";
385  return;
386  }
387 
388  // Get RecoMuons
390  if (!event.getByToken(muonToken_, muonHandle)) {
391  edm::LogInfo("RPCRecHitValid") << "Cannot find muon collection\n";
392  return;
393  }
394 
395  typedef edm::PSimHitContainer::const_iterator SimHitIter;
396  typedef RPCRecHitCollection::const_iterator RecHitIter;
397  typedef std::vector<TrackPSimHitRef> SimHitRefs;
398 
399  // TrackingParticles with (and without) RPC simHits
400  SimHitRefs muonSimHits, pthrSimHits;
401 
402  for (int i = 0, n = simParticleHandle->size(); i < n; ++i) {
403  TrackingParticleRef simParticle(simParticleHandle, i);
404  if (simParticle->pt() < 1.0 or simParticle->p() < 2.5)
405  continue; // globalMuon acceptance
406 
407  // Collect SimHits from this Tracking Particle
408  SimHitRefs simHitsFromParticle;
409  auto range = std::equal_range(simHitsTPAssoc->begin(),
410  simHitsTPAssoc->end(),
411  std::make_pair(simParticle, TrackPSimHitRef()),
413  for (auto simParticleToHit = range.first; simParticleToHit != range.second; ++simParticleToHit) {
414  auto simHit = simParticleToHit->second;
415  const DetId detId(simHit->detUnitId());
416  if (detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC)
417  continue;
418 
419  simHitsFromParticle.push_back(simParticleToHit->second);
420  }
421  const int nRPCHit = simHitsFromParticle.size();
422  const bool hasRPCHit = nRPCHit > 0;
423 
424  if (abs(simParticle->pdgId()) == 13) {
425  muonSimHits.insert(muonSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
426 
427  // Count number of Barrel hits and Endcap hits
428  int nRPCHitBarrel = 0;
429  int nRPCHitEndcap = 0;
430  for (const auto &simHit : simHitsFromParticle) {
431  const RPCDetId rpcDetId = static_cast<const RPCDetId>(simHit->detUnitId());
432  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rpcDetId));
433  if (!roll)
434  continue;
435 
436  if (rpcDetId.region() == 0)
437  ++nRPCHitBarrel;
438  else
439  ++nRPCHitEndcap;
440  }
441 
442  // Fill TrackingParticle related histograms
443  h_nRPCHitPerSimMuon->Fill(nRPCHit);
444  if (nRPCHitBarrel and nRPCHitEndcap) {
445  h_nRPCHitPerSimMuonOverlap->Fill(nRPCHit);
446  h_simMuonOverlap_pt->Fill(simParticle->pt());
447  h_simMuonOverlap_eta->Fill(simParticle->eta());
448  h_simMuonOverlap_phi->Fill(simParticle->phi());
449  } else if (nRPCHitBarrel) {
450  h_nRPCHitPerSimMuonBarrel->Fill(nRPCHit);
451  h_simMuonBarrel_pt->Fill(simParticle->pt());
452  h_simMuonBarrel_eta->Fill(simParticle->eta());
453  h_simMuonBarrel_phi->Fill(simParticle->phi());
454  } else if (nRPCHitEndcap) {
455  h_nRPCHitPerSimMuonEndcap->Fill(nRPCHit);
456  h_simMuonEndcap_pt->Fill(simParticle->pt());
457  h_simMuonEndcap_eta->Fill(simParticle->eta());
458  h_simMuonEndcap_phi->Fill(simParticle->phi());
459  } else {
460  h_simMuonNoRPC_pt->Fill(simParticle->pt());
461  h_simMuonNoRPC_eta->Fill(simParticle->eta());
462  h_simMuonNoRPC_phi->Fill(simParticle->phi());
463  }
464  } else {
465  pthrSimHits.insert(pthrSimHits.end(), simHitsFromParticle.begin(), simHitsFromParticle.end());
466  }
467 
468  if (hasRPCHit) {
469  switch (simParticle->pdgId()) {
470  case 13:
471  h_simParticleType->Fill(0);
472  break;
473  case -13:
474  h_simParticleType->Fill(1);
475  break;
476  case 11:
477  h_simParticleType->Fill(2);
478  break;
479  case -11:
480  h_simParticleType->Fill(3);
481  break;
482  case 211:
483  h_simParticleType->Fill(4);
484  break;
485  case -211:
486  h_simParticleType->Fill(5);
487  break;
488  case 321:
489  h_simParticleType->Fill(6);
490  break;
491  case -321:
492  h_simParticleType->Fill(7);
493  break;
494  case 2212:
495  h_simParticleType->Fill(8);
496  break;
497  case -2212:
498  h_simParticleType->Fill(9);
499  break;
500  default:
501  h_simParticleType->Fill(10);
502  break;
503  }
504  }
505  }
506 
507  // Loop over muon simHits, fill histograms which does not need associations
508  int nRefHitBarrel = 0, nRefHitEndcap = 0;
509  for (const auto &simHit : muonSimHits) {
510  const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
511  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
512 
513  const int region = roll->id().region();
514  const int ring = roll->id().ring();
515  // const int sector = roll->id().sector();
516  const int station = roll->id().station();
517  // const int layer = roll->id().layer();
518  // const int subSector = roll->id().subsector();
519 
520  if (region == 0) {
521  ++nRefHitBarrel;
522  h_.refHitOccupancyBarrel_wheel->Fill(ring);
523  h_.refHitOccupancyBarrel_station->Fill(station);
524  h_.refHitOccupancyBarrel_wheel_station->Fill(ring, station);
525 
526  h_refOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[simHit->detUnitId()]);
527  } else {
528  ++nRefHitEndcap;
529  h_.refHitOccupancyEndcap_disk->Fill(region * station);
530  h_.refHitOccupancyEndcap_disk_ring->Fill(region * station, ring);
531 
532  h_refOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[simHit->detUnitId()]);
533  }
534  }
535 
536  // Loop over punch-through simHits, fill histograms which does not need
537  // associations
538  for (const auto &simHit : pthrSimHits) {
539  const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
540  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
541 
542  const int region = roll->id().region();
543  const int ring = roll->id().ring();
544  // const int sector = roll->id().sector();
545  const int station = roll->id().station();
546  // const int layer = roll->id().layer();
547  // const int subSector = roll->id().subsector();
548 
549  if (region == 0) {
550  ++nRefHitBarrel;
551  h_refPunchOccupancyBarrel_wheel->Fill(ring);
552  h_refPunchOccupancyBarrel_station->Fill(station);
553  h_refPunchOccupancyBarrel_wheel_station->Fill(ring, station);
554 
555  h_refOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[simHit->detUnitId()]);
556  } else {
557  ++nRefHitEndcap;
558  h_refPunchOccupancyEndcap_disk->Fill(region * station);
559  h_refPunchOccupancyEndcap_disk_ring->Fill(region * station, ring);
560 
561  h_refOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[simHit->detUnitId()]);
562  }
563  }
564  h_.nRefHitBarrel->Fill(nRefHitBarrel);
565  h_.nRefHitEndcap->Fill(nRefHitEndcap);
566 
567  // Loop over recHits, fill histograms which does not need associations
568  int sumClusterSizeBarrel = 0, sumClusterSizeEndcap = 0;
569  int nRecHitBarrel = 0, nRecHitEndcap = 0;
570  for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
571  const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
572  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
573  if (!roll)
574  continue;
575 
576  const int region = roll->id().region();
577  const int ring = roll->id().ring();
578  // const int sector = roll->id().sector();
579  const int station = roll->id().station();
580  // const int layer = roll->id().layer();
581  // const int subSector = roll->id().subsector();
582 
583  const double time = recHitIter->timeError() >= 0 ? recHitIter->time() : recHitIter->BunchX() * 25;
584 
585  h_.clusterSize->Fill(recHitIter->clusterSize());
586 
587  if (region == 0) {
588  ++nRecHitBarrel;
589  sumClusterSizeBarrel += recHitIter->clusterSize();
590  h_.clusterSizeBarrel->Fill(recHitIter->clusterSize());
591  h_.recHitOccupancyBarrel_wheel->Fill(ring);
592  h_.recHitOccupancyBarrel_station->Fill(station);
593  h_.recHitOccupancyBarrel_wheel_station->Fill(ring, station);
594 
595  h_.timeBarrel->Fill(time);
596  } else {
597  ++nRecHitEndcap;
598  sumClusterSizeEndcap += recHitIter->clusterSize();
599  h_.clusterSizeEndcap->Fill(recHitIter->clusterSize());
600  h_.recHitOccupancyEndcap_disk->Fill(region * station);
601  h_.recHitOccupancyEndcap_disk_ring->Fill(region * station, ring);
602 
603  h_.timeEndcap->Fill(time);
604  }
605 
606  if (roll->isIRPC()) {
607  h_.timeIRPC->Fill(time);
608  } else {
609  h_.timeCRPC->Fill(time);
610  }
611  }
612  const double nRecHit = nRecHitBarrel + nRecHitEndcap;
613  h_.nRecHitBarrel->Fill(nRecHitBarrel);
614  h_.nRecHitEndcap->Fill(nRecHitEndcap);
615  if (nRecHit > 0) {
616  const int sumClusterSize = sumClusterSizeBarrel + sumClusterSizeEndcap;
617  h_.avgClusterSize->Fill(double(sumClusterSize) / nRecHit);
618 
619  if (nRecHitBarrel > 0) {
620  h_.avgClusterSizeBarrel->Fill(double(sumClusterSizeBarrel) / nRecHitBarrel);
621  }
622  if (nRecHitEndcap > 0) {
623  h_.avgClusterSizeEndcap->Fill(double(sumClusterSizeEndcap) / nRecHitEndcap);
624  }
625  }
626 
627  // Start matching SimHits to RecHits
628  typedef std::map<TrackPSimHitRef, RecHitIter> SimToRecHitMap;
629  SimToRecHitMap simToRecHitMap;
630 
631  for (const auto &simHit : muonSimHits) {
632  const RPCDetId simDetId = static_cast<const RPCDetId>(simHit->detUnitId());
633  // const RPCRoll* simRoll = dynamic_cast<const
634  // RPCRoll*>(rpcGeom->roll(simDetId));
635 
636  const double simX = simHit->localPosition().x();
637 
638  for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
639  const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
640  const RPCRoll *recRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(recDetId));
641  if (!recRoll)
642  continue;
643 
644  if (simDetId != recDetId)
645  continue;
646 
647  const double recX = recHitIter->localPosition().x();
648  const double newDx = fabs(recX - simX);
649 
650  // Associate SimHit to RecHit
651  SimToRecHitMap::const_iterator prevSimToReco = simToRecHitMap.find(simHit);
652  if (prevSimToReco == simToRecHitMap.end()) {
653  simToRecHitMap.insert(std::make_pair(simHit, recHitIter));
654  } else {
655  const double oldDx = fabs(prevSimToReco->second->localPosition().x() - simX);
656 
657  if (newDx < oldDx) {
658  simToRecHitMap[simHit] = recHitIter;
659  }
660  }
661  }
662  }
663 
664  // Now we have simHit-recHit mapping
665  // So we can fill up relavant histograms
666  int nMatchHitBarrel = 0, nMatchHitEndcap = 0;
667  for (const auto &match : simToRecHitMap) {
668  TrackPSimHitRef simHit = match.first;
669  RecHitIter recHitIter = match.second;
670 
671  const RPCDetId detId = static_cast<const RPCDetId>(simHit->detUnitId());
672  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
673 
674  const int region = roll->id().region();
675  const int ring = roll->id().ring();
676  // const int sector = roll->id().sector();
677  const int station = roll->id().station();
678  // const int layer = roll->id().layer();
679  // const int subsector = roll->id().subsector();
680 
681  const double simX = simHit->localPosition().x();
682  const double recX = recHitIter->localPosition().x();
683  const double errX = sqrt(recHitIter->localPositionError().xx());
684  const double dX = recX - simX;
685  const double pull = errX == 0 ? -999 : dX / errX;
686 
687  // const GlobalPoint simPos = roll->toGlobal(simHitIter->localPosition());
688  // const GlobalPoint recPos = roll->toGlobal(recHitIter->localPosition());
689 
690  if (region == 0) {
691  ++nMatchHitBarrel;
692  h_.resBarrel->Fill(dX);
693  h_.pullBarrel->Fill(pull);
694  h_.matchOccupancyBarrel_wheel->Fill(ring);
695  h_.matchOccupancyBarrel_station->Fill(station);
696  h_.matchOccupancyBarrel_wheel_station->Fill(ring, station);
697 
698  h_.res_wheel_res->Fill(ring, dX);
699  h_.res_station_res->Fill(station, dX);
700  h_.pull_wheel_pull->Fill(ring, pull);
701  h_.pull_station_pull->Fill(station, pull);
702 
703  h_matchOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[detId.rawId()]);
704  } else {
705  ++nMatchHitEndcap;
706  h_.resEndcap->Fill(dX);
707  h_.pullEndcap->Fill(pull);
708  h_.matchOccupancyEndcap_disk->Fill(region * station);
709  h_.matchOccupancyEndcap_disk_ring->Fill(region * station, ring);
710 
711  h_.res_disk_res->Fill(region * station, dX);
712  h_.res_ring_res->Fill(ring, dX);
713  h_.pull_disk_pull->Fill(region * station, pull);
714  h_.pull_ring_pull->Fill(ring, pull);
715 
716  h_matchOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[detId.rawId()]);
717  }
718  }
719  h_.nMatchHitBarrel->Fill(nMatchHitBarrel);
720  h_.nMatchHitEndcap->Fill(nMatchHitEndcap);
721 
722  // Reco Muon hits
723  for (reco::MuonCollection::const_iterator muon = muonHandle->begin(); muon != muonHandle->end(); ++muon) {
724  if (!muon->isGlobalMuon())
725  continue;
726 
727  int nRPCHitBarrel = 0;
728  int nRPCHitEndcap = 0;
729 
730  const reco::TrackRef glbTrack = muon->globalTrack();
731  for (trackingRecHit_iterator recHit = glbTrack->recHitsBegin(); recHit != glbTrack->recHitsEnd(); ++recHit) {
732  if (!(*recHit)->isValid())
733  continue;
734  const DetId detId = (*recHit)->geographicalId();
735  if (detId.det() != DetId::Muon or detId.subdetId() != MuonSubdetId::RPC)
736  continue;
737  const RPCDetId rpcDetId = static_cast<const RPCDetId>(detId);
738 
739  if (rpcDetId.region() == 0)
740  ++nRPCHitBarrel;
741  else
742  ++nRPCHitEndcap;
743  }
744 
745  const int nRPCHit = nRPCHitBarrel + nRPCHitEndcap;
746  h_nRPCHitPerRecoMuon->Fill(nRPCHit);
747  if (nRPCHitBarrel and nRPCHitEndcap) {
748  h_nRPCHitPerRecoMuonOverlap->Fill(nRPCHit);
749  h_recoMuonOverlap_pt->Fill(muon->pt());
750  h_recoMuonOverlap_eta->Fill(muon->eta());
751  h_recoMuonOverlap_phi->Fill(muon->phi());
752  } else if (nRPCHitBarrel) {
753  h_nRPCHitPerRecoMuonBarrel->Fill(nRPCHit);
754  h_recoMuonBarrel_pt->Fill(muon->pt());
755  h_recoMuonBarrel_eta->Fill(muon->eta());
756  h_recoMuonBarrel_phi->Fill(muon->phi());
757  } else if (nRPCHitEndcap) {
758  h_nRPCHitPerRecoMuonEndcap->Fill(nRPCHit);
759  h_recoMuonEndcap_pt->Fill(muon->pt());
760  h_recoMuonEndcap_eta->Fill(muon->eta());
761  h_recoMuonEndcap_phi->Fill(muon->phi());
762  } else {
763  h_recoMuonNoRPC_pt->Fill(muon->pt());
764  h_recoMuonNoRPC_eta->Fill(muon->eta());
765  h_recoMuonNoRPC_phi->Fill(muon->phi());
766  }
767  }
768 
769  // Find Non-muon hits
770  for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
771  const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
772  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
773 
774  const int region = roll->id().region();
775  const int ring = roll->id().ring();
776  // const int sector = roll->id().sector();
777  const int station = roll->id().station();
778  // const int layer = roll->id().layer();
779  // const int subsector = roll->id().subsector();
780 
781  bool matched = false;
782  for (const auto &match : simToRecHitMap) {
783  if (recHitIter == match.second) {
784  matched = true;
785  break;
786  }
787  }
788 
789  if (!matched) {
790  int nPunchMatched = 0;
791  // Check if this recHit came from non-muon simHit
792  for (const auto &simHit : pthrSimHits) {
793  const int absSimHitPType = abs(simHit->particleType());
794  if (absSimHitPType == 13)
795  continue;
796 
797  const RPCDetId simDetId = static_cast<const RPCDetId>(simHit->detUnitId());
798  if (simDetId == detId)
799  ++nPunchMatched;
800  }
801 
802  if (nPunchMatched > 0) {
803  if (region == 0) {
804  h_recPunchOccupancyBarrel_wheel->Fill(ring);
805  h_recPunchOccupancyBarrel_station->Fill(station);
806  h_recPunchOccupancyBarrel_wheel_station->Fill(ring, station);
807  } else {
808  h_recPunchOccupancyEndcap_disk->Fill(region * station);
809  h_recPunchOccupancyEndcap_disk_ring->Fill(region * station, ring);
810  }
811  }
812  }
813  }
814 
815  // Find noise recHits : RecHits without SimHit match
816  for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
817  const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
818  const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(recDetId));
819 
820  const int region = roll->id().region();
821  // const int ring = roll->id().ring(); // UNUSED VARIABLE
822  // const int sector = roll->id().sector();
823  // const int station = roll->id().station(); // UNUSED VARIABLE
824  // const int layer = roll->id().layer();
825  // const int subsector = roll->id().subsector();
826 
827  const double recX = recHitIter->localPosition().x();
828  const double recErrX = sqrt(recHitIter->localPositionError().xx());
829 
830  bool matched = false;
831  for (SimHitIter simHitIter = simHitHandle->begin(); simHitIter != simHitHandle->end(); ++simHitIter) {
832  const RPCDetId simDetId = static_cast<const RPCDetId>(simHitIter->detUnitId());
833  const RPCRoll *simRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(simDetId));
834  if (!simRoll)
835  continue;
836 
837  if (simDetId != recDetId)
838  continue;
839 
840  const double simX = simHitIter->localPosition().x();
841  const double dX = fabs(recX - simX);
842 
843  if (dX / recErrX < 5) {
844  matched = true;
845  break;
846  }
847  }
848 
849  if (!matched) {
850  if (region == 0) {
851  h_noiseOccupancyBarrel_detId->Fill(detIdToIndexMapBarrel_[recDetId.rawId()]);
852  } else {
853  h_noiseOccupancyEndcap_detId->Fill(detIdToIndexMapEndcap_[recDetId.rawId()]);
854  }
855  }
856  }
857 
858  h_eventCount->Fill(2);
859 }
860 
const double Pi
virtual int nstrips() const =0
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
RPCRecHitValid(const edm::ParameterSet &pset)
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
static bool simHitTPAssociationListGreater(SimHitTPPair i, SimHitTPPair j)
virtual TH2F * getTH2F() const
std::vector< const GeomDet * > DetContainer
C::const_iterator const_iterator
constant access iterator type
Definition: RangeMap.h:43
virtual float stripLength() const =0
MonitorElement * bookProfile(TString const &name, TString const &title, int nchX, double lowX, double highX, int, double lowY, double highY, char const *option="s", FUNC onbooking=NOOP())
Definition: DQMStore.h:399
T sqrt(T t)
Definition: SSEVec.h:19
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
RPCRecHitValid::MonitorElement * MEP
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Log< level::Info, false > LogInfo
Definition: DetId.h:17
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:212
static constexpr int RPC
Definition: MuonSubdetId.h:13
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.
Definition: RPCDetId.h:53
virtual TH1 * getTH1() const
virtual float pitch() const =0
edm::Ref< edm::PSimHitContainer > TrackPSimHitRef
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
const std::vector< const RPCRoll * > & rolls() const
Return the Rolls.
Definition: RPCChamber.cc:40
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override
Definition: event.py:1
Definition: Run.h:45