RPCMonitorDigi/src/RPCMonitorDigi.cc

0001 #include "DQM/RPCMonitorDigi/interface/RPCMonitorDigi.h"
0002 #include "DQM/RPCMonitorClient/interface/utils.h"
0003 #include "DQM/RPCMonitorClient/interface/RPCNameHelper.h"
0004
0005 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0006 #include "Geometry/RPCGeometry/interface/RPCGeomServ.h"
0007
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include <set>
0013 #include <fmt/format.h>
0014
0015 const std::array<std::string, 3> RPCMonitorDigi::regionNames_ = {{"Endcap-", "Barrel", "Endcap+"}};
0016
0017 RPCMonitorDigi::RPCMonitorDigi(const edm::ParameterSet& pset)
0018     : counter(0),
0019       muonRPCEvents_(nullptr),
0020       NumberOfRecHitMuon_(nullptr),
0021       NumberOfMuon_(nullptr),
0022       numberOfDisks_(0),
0023       numberOfInnerRings_(0) {
0024   useMuonDigis_ = pset.getUntrackedParameter<bool>("UseMuon", true);
0025   useRollInfo_ = pset.getUntrackedParameter<bool>("UseRollInfo", false);
0026
0027   muPtCut_ = pset.getUntrackedParameter<double>("MuonPtCut", 3.0);
0028   muEtaCut_ = pset.getUntrackedParameter<double>("MuonEtaCut", 1.9);
0029
0030   subsystemFolder_ = pset.getUntrackedParameter<std::string>("RPCFolder", "RPC");
0031   globalFolder_ = pset.getUntrackedParameter<std::string>("GlobalFolder", "SummaryHistograms");
0032
0033   //Parametersets for tokens
0034   muonLabel_ = consumes<reco::CandidateView>(pset.getParameter<edm::InputTag>("MuonLabel"));
0035   rpcRecHitLabel_ = consumes<RPCRecHitCollection>(pset.getParameter<edm::InputTag>("RecHitLabel"));
0036   scalersRawToDigiLabel_ = consumes<DcsStatusCollection>(pset.getParameter<edm::InputTag>("ScalersRawToDigiLabel"));
0037
0038   noiseFolder_ = pset.getUntrackedParameter<std::string>("NoiseFolder", "AllHits");
0039   muonFolder_ = pset.getUntrackedParameter<std::string>("MuonFolder", "Muon");
0040
0041   rpcGeomToken_ = esConsumes<edm::Transition::BeginRun>();
0042 }
0043
0044 void RPCMonitorDigi::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& r, edm::EventSetup const& iSetup) {
0045   edm::LogInfo("rpcmonitordigi") << "[RPCMonitorDigi]: Begin Run ";
0046
0047   numberOfInnerRings_ = 4;  // set default value
0048
0049   std::set<int> disk_set;
0050   const auto& rpcGeo = iSetup.getData(rpcGeomToken_);
0051
0052   //loop on geometry to book all MEs
0053   edm::LogInfo("rpcmonitordigi") << "[RPCMonitorDigi]: Booking histograms per roll. ";
0054   for (auto ich : rpcGeo.dets()) {
0055     const RPCChamber* ch = dynamic_cast<const RPCChamber*>(ich);
0056     if (!ch)
0057       continue;
0058     const auto& roles = ch->rolls();
0059
0060     if (useRollInfo_) {
0061       for (auto roll : roles) {
0062         const RPCDetId& rpcId = roll->id();
0063
0064         //get station and inner ring
0065         if (rpcId.region() != 0) {
0066           disk_set.insert(rpcId.station());
0067           numberOfInnerRings_ = std::min(numberOfInnerRings_, rpcId.ring());
0068         }
0069
0070         //booking all histograms
0071         const std::string nameID = RPCNameHelper::rollName(rpcId);
0072         if (useMuonDigis_)
0073           bookRollME(ibooker, rpcId, &rpcGeo, muonFolder_, meMuonCollection[nameID]);
0074         bookRollME(ibooker, rpcId, &rpcGeo, noiseFolder_, meNoiseCollection[nameID]);
0075       }
0076     } else {
0077       const RPCDetId& rpcId = roles[0]->id();  //any roll would do - here I just take the first one
0078       const std::string nameID = RPCNameHelper::chamberName(rpcId);
0079       if (useMuonDigis_)
0080         bookRollME(ibooker, rpcId, &rpcGeo, muonFolder_, meMuonCollection[nameID]);
0081       bookRollME(ibooker, rpcId, &rpcGeo, noiseFolder_, meNoiseCollection[nameID]);
0082       if (rpcId.region() != 0) {
0083         disk_set.insert(rpcId.station());
0084         numberOfInnerRings_ = std::min(numberOfInnerRings_, rpcId.ring());
0085       }
0086     }
0087   }  //end loop on geometry to book all MEs
0088
0089   numberOfDisks_ = disk_set.size();
0090
0091   //Book
0092   this->bookRegionME(ibooker, noiseFolder_, regionNoiseCollection);
0093   this->bookSectorRingME(ibooker, noiseFolder_, sectorRingNoiseCollection);
0094   this->bookWheelDiskME(ibooker, noiseFolder_, wheelDiskNoiseCollection);
0095
0096   ibooker.setCurrentFolder(subsystemFolder_ + "/" + noiseFolder_);
0097
0098   noiseRPCEvents_ = ibooker.book1D("RPCEvents", "RPCEvents", 1, 0.5, 1.5);
0099
0100   if (useMuonDigis_) {
0101     this->bookRegionME(ibooker, muonFolder_, regionMuonCollection);
0102     this->bookSectorRingME(ibooker, muonFolder_, sectorRingMuonCollection);
0103     this->bookWheelDiskME(ibooker, muonFolder_, wheelDiskMuonCollection);
0104
0105     ibooker.setCurrentFolder(subsystemFolder_ + "/" + muonFolder_);
0106
0107     muonRPCEvents_ = ibooker.book1D("RPCEvents", "RPCEvents", 1, 0.5, 1.5);
0108     NumberOfMuon_ = ibooker.book1D("NumberOfMuons", "Number of Muons", 11, -0.5, 10.5);
0109     NumberOfRecHitMuon_ = ibooker.book1D("NumberOfRecHitMuons", "Number of RPC RecHits per Muon", 8, -0.5, 7.5);
0110   }
0111 }
0112
0113 void RPCMonitorDigi::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0114   counter++;
0115   edm::LogInfo("rpcmonitordigi") << "[RPCMonitorDigi]: Beginning analyzing event " << counter;
0116
0117   //Muons
0118   edm::Handle<reco::CandidateView> muonCands;
0119   event.getByToken(muonLabel_, muonCands);
0120
0121   std::map<RPCDetId, std::vector<RPCRecHit> > rechitMuon;
0122
0123   int numMuons = 0;
0124   int numRPCRecHit = 0;
0125
0126   if (muonCands.isValid()) {
0127     int nStaMuons = muonCands->size();
0128
0129     for (int i = 0; i < nStaMuons; i++) {
0130       const reco::Candidate& goodMuon = (*muonCands)[i];
0131       const reco::Muon* muCand = dynamic_cast<const reco::Muon*>(&goodMuon);
0132
0133       if (!muCand->isGlobalMuon())
0134         continue;
0135       if (muCand->pt() < muPtCut_ || fabs(muCand->eta()) > muEtaCut_)
0136         continue;
0137       numMuons++;
0138       reco::Track muTrack = (*(muCand->outerTrack()));
0139       std::vector<TrackingRecHitRef> rpcTrackRecHits;
0140       //loop on mu rechits
0141       for (trackingRecHit_iterator it = muTrack.recHitsBegin(); it != muTrack.recHitsEnd(); it++) {
0142         if (!(*it)->isValid())
0143           continue;
0144         int muSubDetId = (*it)->geographicalId().subdetId();
0145         if (muSubDetId == MuonSubdetId::RPC) {
0146           numRPCRecHit++;
0147           TrackingRecHit* tkRecHit = (*it)->clone();
0148           RPCRecHit* rpcRecHit = dynamic_cast<RPCRecHit*>(tkRecHit);
0149           int detId = (int)rpcRecHit->rpcId();
0150           if (rechitMuon.find(detId) == rechitMuon.end() || rechitMuon[detId].empty()) {
0151             std::vector<RPCRecHit> myVect(1, *rpcRecHit);
0152             rechitMuon[detId] = myVect;
0153           } else {
0154             rechitMuon[detId].push_back(*rpcRecHit);
0155           }
0156         }
0157       }  // end loop on mu rechits
0158     }
0159
0160     //Fill muon counter
0161     if (NumberOfMuon_) {
0162       NumberOfMuon_->Fill(numMuons);
0163     }
0164
0165     //Fill rechit counter for muons
0166     if (NumberOfRecHitMuon_ && numMuons > 0) {
0167       NumberOfRecHitMuon_->Fill(numRPCRecHit);
0168     }
0169
0170     //Fill counter of RPC events with rechits associated in with a muon
0171     if (muonRPCEvents_ != nullptr && numRPCRecHit > 0) {
0172       muonRPCEvents_->Fill(1);
0173     }
0174
0175     //Perform client operation
0176     this->performSourceOperation(rechitMuon, muonFolder_);
0177
0178   } else {
0179     edm::LogError("rpcmonitordigi") << "[RPCMonitorDigi]: Muons - Product not valid for event" << counter;
0180   }
0181
0182   //RecHits
0183   edm::Handle<RPCRecHitCollection> rpcHits;
0184   event.getByToken(rpcRecHitLabel_, rpcHits);
0185   std::map<RPCDetId, std::vector<RPCRecHit> > rechitNoise;
0186
0187   if (rpcHits.isValid()) {
0188     //    RPC rec hits NOT associated to a muon
0189     for (auto rpcRecHitIter = rpcHits->begin(); rpcRecHitIter != rpcHits->end(); rpcRecHitIter++) {
0190       const RPCRecHit& rpcRecHit = (*rpcRecHitIter);
0191       int detId = (int)rpcRecHit.rpcId();
0192       if (rechitNoise.find(detId) == rechitNoise.end() || rechitNoise[detId].empty()) {
0193         std::vector<RPCRecHit> myVect(1, rpcRecHit);
0194         rechitNoise[detId] = myVect;
0195       } else {
0196         rechitNoise[detId].push_back(rpcRecHit);
0197       }
0198     }
0199   } else {
0200     edm::LogError("rpcmonitordigi") << "[RPCMonitorDigi]: RPCRecHits - Product not valid for event" << counter;
0201   }
0202
0203   //Fill counter for all RPC events
0204   if (noiseRPCEvents_ != nullptr && !rechitNoise.empty()) {
0205     noiseRPCEvents_->Fill(1);
0206   }
0207   //Perform client operation
0208   this->performSourceOperation(rechitNoise, noiseFolder_);
0209 }
0210
0211 void RPCMonitorDigi::performSourceOperation(std::map<RPCDetId, std::vector<RPCRecHit> >& recHitMap,
0212                                             std::string recHittype) {
0213   edm::LogInfo("rpcmonitordigi") << "[RPCMonitorDigi]: Performing DQM source operations for ";
0214
0215   if (recHitMap.empty())
0216     return;
0217
0218   std::map<std::string, std::map<std::string, MonitorElement*> > meRollCollection;
0219   std::map<std::string, MonitorElement*> meWheelDisk;
0220   std::map<std::string, MonitorElement*> meRegion;
0221   std::map<std::string, MonitorElement*> meSectorRing;
0222
0223   if (recHittype == muonFolder_) {
0224     meRollCollection = meMuonCollection;
0225     meWheelDisk = wheelDiskMuonCollection;
0226     meRegion = regionMuonCollection;
0227     meSectorRing = sectorRingMuonCollection;
0228   } else if (recHittype == noiseFolder_) {
0229     meRollCollection = meNoiseCollection;
0230     meWheelDisk = wheelDiskNoiseCollection;
0231     meRegion = regionNoiseCollection;
0232     meSectorRing = sectorRingNoiseCollection;
0233   } else {
0234     edm::LogWarning("rpcmonitordigi") << "[RPCMonitorDigi]: RecHit type not valid.";
0235     return;
0236   }
0237
0238   int totalNumberOfRecHits[3] = {0, 0, 0};
0239
0240   //Loop on Rolls
0241   for (std::map<RPCDetId, std::vector<RPCRecHit> >::const_iterator detIdIter = recHitMap.begin();
0242        detIdIter != recHitMap.end();
0243        detIdIter++) {
0244     RPCDetId detId = (*detIdIter).first;
0245     RPCGeomServ geoServ(detId);
0246
0247     //get roll number
0248     rpcdqm::utils rpcUtils;
0249     int nr = rpcUtils.detId2RollNr(detId);
0250
0251     const std::string nameRoll = RPCNameHelper::name(detId, useRollInfo_);
0252
0253     int region = (int)detId.region();
0254     int wheelOrDiskNumber;
0255     std::string wheelOrDiskType;
0256     int ring = 0;
0257     int sector = detId.sector();
0258     int totalRolls = 3;
0259     int roll = detId.roll();
0260     if (region == 0) {
0261       wheelOrDiskType = "Wheel";
0262       wheelOrDiskNumber = (int)detId.ring();
0263       int station = detId.station();
0264
0265       if (station == 1) {
0266         if (detId.layer() == 1) {
0267           totalRolls = 2;  //RB1in
0268         } else {
0269           totalRolls = 2;  //RB1out
0270         }
0271         if (roll == 3)
0272           roll = 2;  // roll=3 is Forward
0273       } else if (station == 2) {
0274         if (detId.layer() == 1) {
0275           //RB2in
0276           if (abs(wheelOrDiskNumber) == 2 && roll == 3) {
0277             roll = 2;  //W -2, +2 RB2in has only 2 rolls
0278             totalRolls = 2;
0279           }
0280         } else {
0281           //RB2out
0282           if (abs(wheelOrDiskNumber) != 2 && roll == 3) {
0283             roll = 2;  //W -1, 0, +1 RB2out has only 2 rolls
0284             totalRolls = 2;
0285           }
0286         }
0287       } else if (station == 3) {
0288         totalRolls = 2;  //RB3
0289         if (roll == 3)
0290           roll = 2;
0291       } else {
0292         totalRolls = 2;  //RB4
0293         if (roll == 3)
0294           roll = 2;
0295       }
0296
0297     } else {
0298       wheelOrDiskType = "Disk";
0299       wheelOrDiskNumber = region * (int)detId.station();
0300       ring = detId.ring();
0301     }
0302
0303     std::vector<RPCRecHit> recHits = (*detIdIter).second;
0304     const int numberOfRecHits = recHits.size();
0305     totalNumberOfRecHits[region + 1] += numberOfRecHits;
0306
0307     std::set<int> bxSet;
0308     int numDigi = 0;
0309
0310     std::map<std::string, MonitorElement*> meMap = meRollCollection[nameRoll];
0311
0312     //Loop on recHits
0313     std::string tmpName;
0314     for (std::vector<RPCRecHit>::const_iterator recHitIter = recHits.begin(); recHitIter != recHits.end();
0315          recHitIter++) {
0316       const RPCRecHit& recHit = (*recHitIter);
0317
0318       int bx = recHit.BunchX();
0319       bxSet.insert(bx);
0320       int clusterSize = (int)recHit.clusterSize();
0321       numDigi += clusterSize;
0322       int firstStrip = recHit.firstClusterStrip();
0323       int lastStrip = clusterSize + firstStrip - 1;
0324
0325       // ###################### Roll Level  #################################
0326
0327       tmpName = "Occupancy_" + nameRoll;
0328       if (meMap[tmpName]) {
0329         for (int s = firstStrip; s <= lastStrip; s++) {
0330           if (useRollInfo_) {
0331             meMap[tmpName]->Fill(s);
0332           } else {
0333             const int nstrips = meMap[tmpName]->getNbinsX() / totalRolls;
0334             meMap[tmpName]->Fill(s + nstrips * (roll - 1));
0335           }
0336         }
0337       }
0338
0339       tmpName = "BXDistribution_" + nameRoll;
0340       if (meMap[tmpName])
0341         meMap[tmpName]->Fill(bx);
0342
0343       // ###################### Sector- Ring Level #################################
0344
0345       tmpName = fmt::format("Occupancy_{}_{}_Sector_{}", wheelOrDiskType, wheelOrDiskNumber, sector);
0346       if (meSectorRing[tmpName]) {
0347         for (int s = firstStrip; s <= lastStrip; s++) {  //Loop on digis
0348           meSectorRing[tmpName]->Fill(s, nr);
0349         }
0350       }
0351
0352       tmpName = fmt::format("ClusterSize_{}_{}_Sector_{}", wheelOrDiskType, wheelOrDiskNumber, sector);
0353       if (meSectorRing[tmpName]) {
0354         if (clusterSize >= meSectorRing[tmpName]->getNbinsX())
0355           meSectorRing[tmpName]->Fill(meSectorRing[tmpName]->getNbinsX(), nr);
0356         else
0357           meSectorRing[tmpName]->Fill(clusterSize, nr);
0358       }
0359
0360       tmpName = fmt::format("Occupancy_{}_{}_Ring_{}", wheelOrDiskType, wheelOrDiskNumber, ring);
0361       if (geoServ.segment() > 0 && geoServ.segment() <= 18) {
0362         tmpName += "_CH01-CH18";
0363       } else if (geoServ.segment() >= 19) {
0364         tmpName += "_CH19-CH36";
0365       }
0366
0367       if (meSectorRing[tmpName]) {
0368         for (int s = firstStrip; s <= lastStrip; s++) {  //Loop on digis
0369           meSectorRing[tmpName]->Fill(s + 32 * (detId.roll() - 1), geoServ.segment());
0370         }
0371       }
0372
0373       tmpName = fmt::format("ClusterSize_{}_{}_Ring_{}", wheelOrDiskType, wheelOrDiskNumber, ring);
0374       if (geoServ.segment() > 0 && geoServ.segment() <= 9) {
0375         tmpName += "_CH01-CH09";
0376       } else if (geoServ.segment() >= 10 && geoServ.segment() <= 18) {
0377         tmpName += "_CH10-CH18";
0378       } else if (geoServ.segment() >= 19 && geoServ.segment() <= 27) {
0379         tmpName += "_CH19-CH27";
0380       } else if (geoServ.segment() >= 28 && geoServ.segment() <= 36) {
0381         tmpName += "_CH28-CH36";
0382       }
0383
0384       if (meSectorRing[tmpName]) {
0385         if (clusterSize >= meSectorRing[tmpName]->getNbinsX())
0386           meSectorRing[tmpName]->Fill(meSectorRing[tmpName]->getNbinsX(),
0387                                       3 * (geoServ.segment() - 1) + (3 - detId.roll()) + 1);
0388         else
0389           meSectorRing[tmpName]->Fill(clusterSize, 3 * (geoServ.segment() - 1) + (3 - detId.roll()) + 1);
0390       }
0391
0392       // ###################### Wheel/Disk Level #########################
0393       if (region == 0) {
0394         tmpName = fmt::format("1DOccupancy_Wheel_{}", wheelOrDiskNumber);
0395         if (meWheelDisk[tmpName])
0396           meWheelDisk[tmpName]->Fill(sector, clusterSize);
0397
0398         tmpName = fmt::format("Occupancy_Roll_vs_Sector_{}_{}", wheelOrDiskType, wheelOrDiskNumber);
0399         if (meWheelDisk[tmpName])
0400           meWheelDisk[tmpName]->Fill(sector, nr, clusterSize);
0401
0402       } else {
0403         tmpName = fmt::format("1DOccupancy_Ring_{}", ring);
0404         if ((meWheelDisk[tmpName])) {
0405           if (wheelOrDiskNumber > 0) {
0406             meWheelDisk[tmpName]->Fill(wheelOrDiskNumber + numberOfDisks_, clusterSize);
0407           } else {
0408             meWheelDisk[tmpName]->Fill(wheelOrDiskNumber + numberOfDisks_ + 1, clusterSize);
0409           }
0410         }
0411
0412         tmpName = fmt::format("Occupancy_Ring_vs_Segment_{}_{}", wheelOrDiskType, wheelOrDiskNumber);
0413         if (meWheelDisk[tmpName])
0414           meWheelDisk[tmpName]->Fill(geoServ.segment(), (ring - 1) * 3 - detId.roll() + 1, clusterSize);
0415       }
0416
0417       tmpName = fmt::format("BxDistribution_{}_{}", wheelOrDiskType, wheelOrDiskNumber);
0418       if (meWheelDisk[tmpName])
0419         meWheelDisk[tmpName]->Fill(bx);
0420
0421     }  //end loop on recHits
0422
0423     tmpName = "BXWithData_" + nameRoll;
0424     if (meMap[tmpName])
0425       meMap[tmpName]->Fill(bxSet.size());
0426
0427     tmpName = "NumberOfClusters_" + nameRoll;
0428     if (meMap[tmpName])
0429       meMap[tmpName]->Fill(numberOfRecHits);
0430
0431     tmpName = "Multiplicity_" + RPCMonitorDigi::regionNames_[region + 1];
0432     if (meRegion[tmpName])
0433       meRegion[tmpName]->Fill(numDigi);
0434
0435     if (region == 0) {
0436       if (meRegion["Occupancy_for_Barrel"])
0437         meRegion["Occupancy_for_Barrel"]->Fill(sector, wheelOrDiskNumber, numDigi);
0438     } else {
0439       int xbin = wheelOrDiskNumber + numberOfDisks_;
0440       if (region == -1) {
0441         xbin = wheelOrDiskNumber + numberOfDisks_ + 1;
0442       }
0443       if (meRegion["Occupancy_for_Endcap"]) {
0444         meRegion["Occupancy_for_Endcap"]->Fill(xbin, ring, numDigi);
0445       }
0446     }
0447
0448     tmpName = "Multiplicity_" + nameRoll;
0449     if (meMap[tmpName])
0450       meMap[tmpName]->Fill(numDigi);
0451
0452   }  //end loop on rolls
0453
0454   for (int i = 0; i < 3; i++) {
0455     const std::string tmpName = "NumberOfClusters_" + RPCMonitorDigi::regionNames_[i];
0456     if (meRegion[tmpName])
0457       meRegion[tmpName]->Fill(totalNumberOfRecHits[i]);
0458   }
0459 }