RPCRecHits/src/RPCPointVsRecHit.cc

0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "Validation/RPCRecHits/interface/RPCPointVsRecHit.h"
0003
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006
0007 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0008
0009 using namespace std;
0010
0011 typedef RPCPointVsRecHit::MonitorElement *MEP;
0012
0013 RPCPointVsRecHit::RPCPointVsRecHit(const edm::ParameterSet &pset) {
0014   refHitToken_ = consumes<RPCRecHitCollection>(pset.getParameter<edm::InputTag>("refHit"));
0015   recHitToken_ = consumes<RPCRecHitCollection>(pset.getParameter<edm::InputTag>("recHit"));
0016
0017   subDir_ = pset.getParameter<std::string>("subDir");
0018
0019   rpcGeomToken_ = esConsumes();
0020 }
0021
0022 void RPCPointVsRecHit::analyze(const edm::Event &event, const edm::EventSetup &eventSetup) {
0023   // Get the RPC Geometry
0024   auto rpcGeom = eventSetup.getHandle(rpcGeomToken_);
0025
0026   // Retrieve RefHits from the event
0027   edm::Handle<RPCRecHitCollection> refHitHandle;
0028   if (!event.getByToken(refHitToken_, refHitHandle)) {
0029     edm::LogInfo("RPCPointVsRecHit") << "Cannot find reference hit collection\n";
0030     return;
0031   }
0032
0033   // Retrieve RecHits from the event
0034   edm::Handle<RPCRecHitCollection> recHitHandle;
0035   if (!event.getByToken(recHitToken_, recHitHandle)) {
0036     edm::LogInfo("RPCPointVsRecHit") << "Cannot find recHit collection\n";
0037     return;
0038   }
0039
0040   typedef RPCRecHitCollection::const_iterator RecHitIter;
0041
0042   // Loop over refHits, fill histograms which does not need associations
0043   int nRefHitBarrel = 0, nRefHitEndcap = 0;
0044   for (RecHitIter refHitIter = refHitHandle->begin(); refHitIter != refHitHandle->end(); ++refHitIter) {
0045     const RPCDetId detId = static_cast<const RPCDetId>(refHitIter->rpcId());
0046     const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
0047     if (!roll)
0048       continue;
0049
0050     const int region = roll->id().region();
0051     const int ring = roll->id().ring();
0052     // const int sector = roll->id().sector();
0053     const int station = roll->id().station();
0054     // const int layer = roll->id().layer();
0055     // const int subSector = roll->id().subsector();
0056
0057     if (region == 0) {
0058       h_.refHitOccupancyBarrel_wheel->Fill(ring);
0059       h_.refHitOccupancyBarrel_station->Fill(station);
0060       h_.refHitOccupancyBarrel_wheel_station->Fill(ring, station);
0061     } else {
0062       h_.refHitOccupancyEndcap_disk->Fill(region * station);
0063       h_.refHitOccupancyEndcap_disk_ring->Fill(region * station, ring);
0064     }
0065   }
0066   h_.nRefHitBarrel->Fill(nRefHitBarrel);
0067   h_.nRefHitEndcap->Fill(nRefHitEndcap);
0068
0069   // Loop over recHits, fill histograms which does not need associations
0070   int sumClusterSizeBarrel = 0, sumClusterSizeEndcap = 0;
0071   int nRecHitBarrel = 0, nRecHitEndcap = 0;
0072   for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0073     const RPCDetId detId = static_cast<const RPCDetId>(recHitIter->rpcId());
0074     const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId()));
0075     if (!roll)
0076       continue;
0077
0078     const int region = roll->id().region();
0079     const int ring = roll->id().ring();
0080     // const int sector = roll->id().sector();
0081     const int station = (roll->id().station());
0082     // const int layer = roll->id().layer();
0083     // const int subSector = roll->id().subsector();
0084
0085     h_.clusterSize->Fill(recHitIter->clusterSize());
0086
0087     if (region == 0) {
0088       ++nRecHitBarrel;
0089       sumClusterSizeBarrel += recHitIter->clusterSize();
0090       h_.clusterSizeBarrel->Fill(recHitIter->clusterSize());
0091       h_.recHitOccupancyBarrel_wheel->Fill(ring);
0092       h_.recHitOccupancyBarrel_station->Fill(station);
0093       h_.recHitOccupancyBarrel_wheel_station->Fill(ring, station);
0094     } else {
0095       ++nRecHitEndcap;
0096       sumClusterSizeEndcap += recHitIter->clusterSize();
0097       h_.clusterSizeEndcap->Fill(recHitIter->clusterSize());
0098       h_.recHitOccupancyEndcap_disk->Fill(region * station);
0099       h_.recHitOccupancyEndcap_disk_ring->Fill(ring, region * station);
0100     }
0101   }
0102   const double nRecHit = nRecHitBarrel + nRecHitEndcap;
0103   h_.nRecHitBarrel->Fill(nRecHitBarrel);
0104   h_.nRecHitEndcap->Fill(nRecHitEndcap);
0105   if (nRecHit > 0) {
0106     const int sumClusterSize = sumClusterSizeBarrel + sumClusterSizeEndcap;
0107     h_.avgClusterSize->Fill(double(sumClusterSize) / nRecHit);
0108
0109     if (nRecHitBarrel > 0) {
0110       h_.avgClusterSizeBarrel->Fill(double(sumClusterSizeBarrel) / nRecHitBarrel);
0111     }
0112     if (nRecHitEndcap > 0) {
0113       h_.avgClusterSizeEndcap->Fill(double(sumClusterSizeEndcap) / nRecHitEndcap);
0114     }
0115   }
0116
0117   // Start matching RefHits to RecHits
0118   typedef std::map<RecHitIter, RecHitIter> RecToRecHitMap;
0119   RecToRecHitMap refToRecHitMap;
0120
0121   for (RecHitIter refHitIter = refHitHandle->begin(); refHitIter != refHitHandle->end(); ++refHitIter) {
0122     const RPCDetId refDetId = static_cast<const RPCDetId>(refHitIter->rpcId());
0123     const RPCRoll *refRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(refDetId));
0124     if (!refRoll)
0125       continue;
0126
0127     const double refX = refHitIter->localPosition().x();
0128
0129     for (RecHitIter recHitIter = recHitHandle->begin(); recHitIter != recHitHandle->end(); ++recHitIter) {
0130       const RPCDetId recDetId = static_cast<const RPCDetId>(recHitIter->rpcId());
0131       const RPCRoll *recRoll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(recDetId));
0132       if (!recRoll)
0133         continue;
0134
0135       if (refDetId != recDetId)
0136         continue;
0137
0138       const double recX = recHitIter->localPosition().x();
0139       const double newDx = fabs(recX - refX);
0140
0141       // Associate RefHit to RecHit
0142       RecToRecHitMap::const_iterator prevRefToReco = refToRecHitMap.find(refHitIter);
0143       if (prevRefToReco == refToRecHitMap.end()) {
0144         refToRecHitMap.insert(std::make_pair(refHitIter, recHitIter));
0145       } else {
0146         const double oldDx = fabs(prevRefToReco->second->localPosition().x() - refX);
0147
0148         if (newDx < oldDx) {
0149           refToRecHitMap[refHitIter] = recHitIter;
0150         }
0151       }
0152     }
0153   }
0154
0155   // Now we have refHit-recHit mapping
0156   // So we can fill up relavant histograms
0157   for (RecToRecHitMap::const_iterator match = refToRecHitMap.begin(); match != refToRecHitMap.end(); ++match) {
0158     RecHitIter refHitIter = match->first;
0159     RecHitIter recHitIter = match->second;
0160
0161     const RPCDetId detId = static_cast<const RPCDetId>(refHitIter->rpcId());
0162     const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(detId));
0163
0164     const int region = roll->id().region();
0165     const int ring = roll->id().ring();
0166     // const int sector = roll->id().sector();
0167     const int station = roll->id().station();
0168     // const int layer = roll->id().layer();
0169     // const int subsector = roll->id().subsector();
0170
0171     const double refX = refHitIter->localPosition().x();
0172     const double recX = recHitIter->localPosition().x();
0173     const double errX = sqrt(recHitIter->localPositionError().xx());
0174     const double dX = recX - refX;
0175     const double pull = errX == 0 ? -999 : dX / errX;
0176
0177     // const GlobalPoint refPos = roll->toGlobal(refHitIter->localPosition());
0178     // const GlobalPoint recPos = roll->toGlobal(recHitIter->localPosition());
0179
0180     if (region == 0) {
0181       h_.resBarrel->Fill(dX);
0182       h_.pullBarrel->Fill(pull);
0183       h_.matchOccupancyBarrel_wheel->Fill(ring);
0184       h_.matchOccupancyBarrel_station->Fill(station);
0185       h_.matchOccupancyBarrel_wheel_station->Fill(ring, station);
0186
0187       h_.res_wheel_res->Fill(ring, dX);
0188       h_.res_station_res->Fill(station, dX);
0189       h_.pull_wheel_pull->Fill(ring, pull);
0190       h_.pull_station_pull->Fill(station, pull);
0191     } else {
0192       h_.resEndcap->Fill(dX);
0193       h_.pullEndcap->Fill(pull);
0194       h_.matchOccupancyEndcap_disk->Fill(region * station);
0195       h_.matchOccupancyEndcap_disk_ring->Fill(region * station, ring);
0196
0197       h_.res_disk_res->Fill(region * station, dX);
0198       h_.res_ring_res->Fill(ring, dX);
0199       h_.pull_disk_pull->Fill(region * station, pull);
0200       h_.pull_ring_pull->Fill(ring, pull);
0201     }
0202   }
0203 }
0204
0205 void RPCPointVsRecHit::bookHistograms(DQMStore::IBooker &booker,
0206                                       edm::Run const &run,
0207                                       edm::EventSetup const &eventSetup) {
0208   // Book MonitorElements
0209   h_.bookHistograms(booker, subDir_);
0210 }
0211
0212 DEFINE_FWK_MODULE(RPCPointVsRecHit);