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