MuonRPCDigis/src/RPCDigiValid.cc

0001 #include "Validation/MuonRPCDigis/interface/RPCDigiValid.h"
0002
0003 #include "FWCore/Utilities/interface/InputTag.h"
0004 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0005
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0008 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0009 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0010 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0011
0012 #include <cmath>
0013 #include <fmt/format.h>
0014
0015 using namespace std;
0016 using namespace edm;
0017
0018 RPCDigiValid::RPCDigiValid(const ParameterSet &ps) {
0019   //  Init the tokens for run data retrieval - stanislav
0020   //  ps.getUntackedParameter<InputTag> retrieves a InputTag from the
0021   //  configuration. The second param is default value module, instance and
0022   //  process labels may be passed in a single string if separated by colon ':'
0023   //  (@see the edm::InputTag constructor documentation)
0024   simHitToken_ = consumes<PSimHitContainer>(
0025       ps.getUntrackedParameter<edm::InputTag>("simHitTag", edm::InputTag("g4SimHits:MuonRPCHits")));
0026   rpcDigiToken_ = consumes<RPCDigiCollection>(
0027       ps.getUntrackedParameter<edm::InputTag>("rpcDigiTag", edm::InputTag("simMuonRPCDigis")));
0028
0029   isDigiTimeAvailable_ = ps.getUntrackedParameter<bool>("digiTime", false);
0030
0031   rpcGeomToken_ = esConsumes();
0032 }
0033
0034 void RPCDigiValid::analyze(const Event &event, const EventSetup &eventSetup) {
0035   // Get the RPC Geometry
0036   auto rpcGeom = eventSetup.getHandle(rpcGeomToken_);
0037
0038   edm::Handle<PSimHitContainer> simHitHandle;
0039   edm::Handle<RPCDigiCollection> rpcDigisHandle;
0040   event.getByToken(simHitToken_, simHitHandle);
0041   event.getByToken(rpcDigiToken_, rpcDigisHandle);
0042
0043   // loop over Simhit
0044   std::map<const RPCRoll *, std::vector<double>> detToSimHitXsMap;
0045   for (auto simIt = simHitHandle->begin(); simIt != simHitHandle->end(); ++simIt) {
0046     const RPCDetId rsid = simIt->detUnitId();
0047     const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rsid));
0048     if (!roll)
0049       continue;
0050
0051     if (detToSimHitXsMap.find(roll) == detToSimHitXsMap.end())
0052       detToSimHitXsMap[roll] = std::vector<double>();
0053     detToSimHitXsMap[roll].push_back(simIt->localPosition().x());
0054
0055     const int region = rsid.region();
0056     const GlobalPoint gp = roll->toGlobal(simIt->localPosition());
0057
0058     hRZ_->Fill(gp.z(), gp.perp());
0059
0060     if (region == 0) {
0061       // Barrel
0062       hXY_Barrel_->Fill(gp.x(), gp.y());
0063
0064       const int station = rsid.station();
0065       const int layer = rsid.layer();
0066       const int stla = (station <= 2) ? (2 * (station - 1) + layer) : (station + 2);
0067
0068       auto match = hZPhi_.find(stla);
0069       if (match != hZPhi_.end()) {
0070         const double phiInDeg = 180. * gp.barePhi() / TMath::Pi();
0071         match->second->Fill(gp.z(), phiInDeg);
0072       }
0073     } else {
0074       // Endcap
0075       const int disk = region * rsid.station();
0076       auto match = hXY_Endcap_.find(disk);
0077       if (match != hXY_Endcap_.end())
0078         match->second->Fill(gp.x(), gp.y());
0079     }
0080   }
0081   for (const auto &detToSimHitXs : detToSimHitXsMap) {
0082     hNSimHitPerRoll_->Fill(detToSimHitXs.second.size());
0083   }
0084
0085   // loop over Digis
0086   std::map<const RPCRoll *, std::vector<double>> detToDigiXsMap;
0087   for (auto detUnitIt = rpcDigisHandle->begin(); detUnitIt != rpcDigisHandle->end(); ++detUnitIt) {
0088     const RPCDetId rsid = (*detUnitIt).first;
0089     const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(rsid));
0090     if (!roll)
0091       continue;
0092     const int region = rsid.region();
0093
0094     const RPCDigiCollection::Range &range = (*detUnitIt).second;
0095
0096     for (auto digiIt = range.first; digiIt != range.second; ++digiIt) {
0097       // Strip profile
0098       const int strip = digiIt->strip();
0099       hStripProf_->Fill(strip);
0100
0101       if (region == 0) {
0102         const int station = rsid.station();
0103         if (station == 1 or station == 2)
0104           hStripProf_RB12_->Fill(strip);
0105         else if (station == 3 or station == 4)
0106           hStripProf_RB34_->Fill(strip);
0107       } else {
0108         if (roll->isIRPC())
0109           hStripProf_IRPC_->Fill(strip);
0110         else
0111           hStripProf_Endcap_->Fill(strip);
0112       }
0113
0114       // Bunch crossing
0115       const int bx = digiIt->bx();
0116       hBxDist_->Fill(bx);
0117       // bx for 4 endcaps
0118       if (rsid.station() == 4) {
0119         if (region == 1) {
0120           hBxDisc_4Plus_->Fill(bx);
0121         } else if (region == -1) {
0122           hBxDisc_4Min_->Fill(bx);
0123         }
0124       }
0125
0126       // Fill timing information
0127       if (isDigiTimeAvailable_) {
0128         const double digiTime = digiIt->hasTime() ? digiIt->time() : digiIt->bx() * 25;
0129         hDigiTimeAll_->Fill(digiTime);
0130         if (digiIt->hasTime()) {
0131           hDigiTime_->Fill(digiTime);
0132           if (roll->isIRPC())
0133             hDigiTimeIRPC_->Fill(digiTime);
0134           else
0135             hDigiTimeNoIRPC_->Fill(digiTime);
0136         }
0137       }
0138
0139       // Keep digi position
0140       const double digiX = roll->centreOfStrip(digiIt->strip()).x();
0141       if (detToDigiXsMap.find(roll) == detToDigiXsMap.end())
0142         detToDigiXsMap[roll] = std::vector<double>();
0143       detToDigiXsMap[roll].push_back(digiX);
0144     }
0145   }
0146   for (const auto &detToDigiXs : detToDigiXsMap) {
0147     const auto digiXs = detToDigiXs.second;
0148     const int nDigi = digiXs.size();
0149     hNDigiPerRoll_->Fill(nDigi);
0150
0151     // Fill residual plots, only for nDigi==1 and nSimHit==1
0152     const auto roll = detToDigiXs.first;
0153     const auto detId = roll->id();
0154     if (nDigi != 1)
0155       continue;
0156     if (detToSimHitXsMap.find(roll) == detToSimHitXsMap.end())
0157       continue;
0158
0159     const auto simHitXs = detToSimHitXsMap[roll];
0160     const int nSimHit = simHitXs.size();
0161     if (nSimHit != 1)
0162       continue;
0163
0164     const double dx = digiXs[0] - simHitXs[0];
0165     hRes_->Fill(dx);
0166     if (roll->isBarrel()) {
0167       const int wheel = detId.ring();  // ring() is wheel number for Barrel
0168       const int station = detId.station();
0169       const int layer = detId.layer();
0170       const int stla = (station <= 2) ? (2 * (station - 1) + layer) : (station + 2);
0171
0172       auto matchLayer = hResBarrelLayers_.find(stla);
0173       if (matchLayer != hResBarrelLayers_.end())
0174         matchLayer->second->Fill(dx);
0175
0176       auto matchWheel = hResBarrelWheels_.find(wheel);
0177       if (matchWheel != hResBarrelWheels_.end())
0178         matchWheel->second->Fill(dx);
0179     } else {
0180       const int disk = detId.region() * detId.station();
0181       auto matchDisk = hResEndcapDisks_.find(disk);
0182       if (matchDisk != hResEndcapDisks_.end())
0183         matchDisk->second->Fill(dx);
0184
0185       auto matchRing = hResEndcapRings_.find(detId.ring());
0186       if (matchRing != hResEndcapRings_.end())
0187         matchRing->second->Fill(dx);
0188     }
0189   }
0190 }
0191
0192 void RPCDigiValid::bookHistograms(DQMStore::IBooker &booker, edm::Run const &run, edm::EventSetup const &eSetup) {
0193   booker.setCurrentFolder("RPCDigisV/RPCDigis");
0194
0195   // Define binnings of 2D-histograms
0196   const double maxZ = 1100;
0197   const int nbinsZ = 220;  // bin width: 10cm
0198   const double maxXY = 800;
0199   const int nbinsXY = 160;  // bin width: 10cm
0200   const double minR = 100, maxR = 800;
0201   const int nbinsR = 70;    // bin width: 10cm
0202   const int nbinsPhi = 72;  // bin width: 5 degree
0203   const double maxBarrelZ = 700;
0204   const int nbinsBarrelZ = 140;  // bin width: 10cm
0205
0206   // RZ plot
0207   hRZ_ = booker.book2D("RZ", "R-Z view;Z (cm);R (cm)", nbinsZ, -maxZ, maxZ, nbinsR, minR, maxR);
0208   hRZ_->setOption("colz");
0209
0210   // XY plots
0211   hXY_Barrel_ = booker.book2D("XY_Barrel", "X-Y view of Barrel", nbinsXY, -maxXY, maxXY, nbinsXY, -maxXY, maxXY);
0212   hXY_Barrel_->setOption("colz");
0213   for (int disk = 1; disk <= 4; ++disk) {
0214     const std::string meNameP = fmt::format("XY_Endcap_p{:1d}", disk);
0215     const std::string meNameN = fmt::format("XY_Endcap_m{:1d}", disk);
0216     const std::string meTitleP = fmt::format("X-Y view of Endcap{:+1d};X (cm);Y (cm)", disk);
0217     const std::string meTitleN = fmt::format("X-Y view of Endcap{:+1d};X (cm);Y (cm)", -disk);
0218     hXY_Endcap_[disk] = booker.book2D(meNameP, meTitleP, nbinsXY, -maxXY, maxXY, nbinsXY, -maxXY, maxXY);
0219     hXY_Endcap_[-disk] = booker.book2D(meNameN, meTitleN, nbinsXY, -maxXY, maxXY, nbinsXY, -maxXY, maxXY);
0220     hXY_Endcap_[disk]->setOption("colz");
0221     hXY_Endcap_[-disk]->setOption("colz");
0222   }
0223
0224   // Z-phi plots
0225   for (int layer = 1; layer <= 6; ++layer) {
0226     const std::string meName = fmt::format("ZPhi_Layer{:1d}", layer);
0227     const std::string meTitle = fmt::format("Z-#phi view of Layer{:1d};Z (cm);#phi (degree)", layer);
0228     hZPhi_[layer] = booker.book2D(meName, meTitle, nbinsBarrelZ, -maxBarrelZ, maxBarrelZ, nbinsPhi, -180, 180);
0229     hZPhi_[layer]->setOption("colz");
0230   }
0231
0232   // Strip profile
0233   hStripProf_ = booker.book1D("Strip_Profile", "Strip_Profile;Strip Number", 100, 0, 100);
0234   hStripProf_RB12_ = booker.book1D("Strip_Profile_RB12", "Strip Profile RB1 and RB2;Strip Number", 92, 0, 92);
0235   hStripProf_RB34_ = booker.book1D("Strip_Profile_RB34", "Strip Profile RB3 and RB4;Strip Number", 62, 0, 62);
0236   hStripProf_Endcap_ = booker.book1D("Strip_Profile_Endcap", "Strip Profile Endcap;Strip Number", 40, 0, 40);
0237   hStripProf_IRPC_ = booker.book1D("Strip_Profile_IRPC", "Strip Profile IRPC;Strip Number", 100, 0, 100);
0238
0239   // Bunch crossing
0240   hBxDist_ = booker.book1D("Bunch_Crossing", "Bunch Crossing;Bunch crossing", 20, -10., 10.);
0241   hBxDisc_4Plus_ = booker.book1D("BxDisc_4Plus", "BxDisc_4Plus", 20, -10., 10.);
0242   hBxDisc_4Min_ = booker.book1D("BxDisc_4Min", "BxDisc_4Min", 20, -10., 10.);
0243
0244   // Timing informations
0245   if (isDigiTimeAvailable_) {
0246     hDigiTimeAll_ =
0247         booker.book1D("DigiTimeAll", "Digi time including present electronics;Digi time (ns)", 100, -12.5, 12.5);
0248     hDigiTime_ = booker.book1D("DigiTime", "Digi time only with timing information;Digi time (ns)", 100, -12.5, 12.5);
0249     hDigiTimeIRPC_ = booker.book1D("DigiTimeIRPC", "IRPC Digi time;Digi time (ns)", 100, -12.5, 12.5);
0250     hDigiTimeNoIRPC_ = booker.book1D("DigiTimeNoIRPC", "non-IRPC Digi time;Digi time (ns)", 100, -12.5, 12.5);
0251   }
0252
0253   // SimHit and Digi multiplicity per roll
0254   hNSimHitPerRoll_ = booker.book1D("NSimHitPerRoll", "SimHit multiplicity per Roll;Multiplicity", 10, 0, 10);
0255   hNDigiPerRoll_ = booker.book1D("NDigiPerRoll", "Digi multiplicity per Roll;Multiplicity", 10, 0, 10);
0256
0257   // Residual of SimHit-Digi x-position
0258   hRes_ = booker.book1D("Digi_SimHit_Difference", "Digi-SimHit difference;dx (cm)", 100, -8, 8);
0259
0260   for (int layer = 1; layer <= 6; ++layer) {
0261     const std::string meName = fmt::format("Residual_Barrel_Layer{:1d}", layer);
0262     const std::string meTitle = fmt::format("Residual of Barrel Layer{:1d};dx (cm)", layer);
0263     hResBarrelLayers_[layer] = booker.book1D(meName, meTitle, 100, -8, 8);
0264   }
0265
0266   hResBarrelWheels_[-2] = booker.book1D("Residual_Barrel_Wheel_m2", "Residual of Barrel Wheel-2;dx (cm)", 100, -8, 8);
0267   hResBarrelWheels_[-1] = booker.book1D("Residual_Barrel_Wheel_m1", "Residual of Barrel Wheel-1;dx (cm)", 100, -8, 8);
0268   hResBarrelWheels_[+0] = booker.book1D("Residual_Barrel_Wheel_00", "Residual of Barrel Wheel 0;dx (cm)", 100, -8, 8);
0269   hResBarrelWheels_[+1] = booker.book1D("Residual_Barrel_Wheel_p1", "Residual of Barrel Wheel+1;dx (cm)", 100, -8, 8);
0270   hResBarrelWheels_[+2] = booker.book1D("Residual_Barrel_Wheel_p2", "Residual of Barrel Wheel+2;dx (cm)", 100, -8, 8);
0271
0272   for (int disk = 1; disk <= 4; ++disk) {
0273     const std::string meNameP = fmt::format("Residual_Endcap_Disk_p{:1d}", disk);
0274     const std::string meNameN = fmt::format("Residual_Endcap_Disk_m{:1d}", disk);
0275     const std::string meTitleP = fmt::format("Residual of Endcap Disk{:+1d};dx (cm)", disk);
0276     const std::string meTitleN = fmt::format("Residual of Endcap Disk{:+1d};dx (cm)", -disk);
0277     hResEndcapDisks_[+disk] = booker.book1D(meNameP, meTitleP, 100, -8, 8);
0278     hResEndcapDisks_[-disk] = booker.book1D(meNameN, meTitleN, 100, -8, 8);
0279   }
0280
0281   hResEndcapRings_[1] = booker.book1D("Residual_Endcap_Ring1", "Residual of Endcap Ring1;dx (cm)", 100, -12, 12);
0282   hResEndcapRings_[2] = booker.book1D("Residual_Endcap_Ring2", "Residual of Endcap Ring2;dx (cm)", 100, -8, 8);
0283   hResEndcapRings_[3] = booker.book1D("Residual_Endcap_Ring3", "Residual of Endcap Ring3;dx (cm)", 100, -8, 8);
0284 }