CMS 3D CMS Logo

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