MuonRPCDigis/src/RPCDigiValid.cc

0001 #include "FWCore/Utilities/interface/InputTag.h"
0002 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0003 #include "Validation/MuonRPCDigis/interface/RPCDigiValid.h"
0004
0005 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0006 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0007 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0008 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0009 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0010
0011 #include <cmath>
0012
0013 using namespace std;
0014 using namespace edm;
0015
0016 RPCDigiValid::RPCDigiValid(const ParameterSet &ps) {
0017   //  Init the tokens for run data retrieval - stanislav
0018   //  ps.getUntackedParameter<InputTag> retrieves a InputTag from the
0019   //  configuration. The second param is default value module, instance and
0020   //  process labels may be passed in a single string if separated by colon ':'
0021   //  (@see the edm::InputTag constructor documentation)
0022   simHitToken = consumes<PSimHitContainer>(
0023       ps.getUntrackedParameter<edm::InputTag>("simHitTag", edm::InputTag("g4SimHits:MuonRPCHits")));
0024   rpcDigiToken = consumes<RPCDigiCollection>(
0025       ps.getUntrackedParameter<edm::InputTag>("rpcDigiTag", edm::InputTag("simMuonRPCDigis")));
0026
0027   outputFile_ = ps.getUntrackedParameter<string>("outputFile", "rpcDigiValidPlots.root");
0028
0029   rpcGeomToken_ = esConsumes();
0030 }
0031
0032 RPCDigiValid::~RPCDigiValid() {}
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> simHit;
0039   edm::Handle<RPCDigiCollection> rpcDigis;
0040   event.getByToken(simHitToken, simHit);
0041   event.getByToken(rpcDigiToken, rpcDigis);
0042
0043   // Loop on simhits
0044   PSimHitContainer::const_iterator simIt;
0045
0046   // loop over Simhit
0047   std::map<RPCDetId, std::vector<double>> allsims;
0048
0049   for (simIt = simHit->begin(); simIt != simHit->end(); simIt++) {
0050     RPCDetId Rsid = (RPCDetId)(*simIt).detUnitId();
0051     const RPCRoll *soll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(Rsid));
0052     int ptype = simIt->particleType();
0053
0054     if (ptype == 13 || ptype == -13) {
0055       std::vector<double> buff;
0056       if (allsims.find(Rsid) != allsims.end()) {
0057         buff = allsims[Rsid];
0058       }
0059
0060       buff.push_back(simIt->localPosition().x());
0061
0062       allsims[Rsid] = buff;
0063     }
0064     GlobalPoint p = soll->toGlobal(simIt->localPosition());
0065
0066     double sim_x = p.x();
0067     double sim_y = p.y();
0068
0069     xyview->Fill(sim_x, sim_y);
0070
0071     if (Rsid.region() == (+1)) {
0072       if (Rsid.station() == 4) {
0073         xyvDplu4->Fill(sim_x, sim_y);
0074       }
0075     } else if (Rsid.region() == (-1)) {
0076       if (Rsid.station() == 4) {
0077         xyvDmin4->Fill(sim_x, sim_y);
0078       }
0079     }
0080     rzview->Fill(p.z(), p.perp());
0081   }
0082   // loop over Digis
0083   RPCDigiCollection::DigiRangeIterator detUnitIt;
0084   for (detUnitIt = rpcDigis->begin(); detUnitIt != rpcDigis->end(); ++detUnitIt) {
0085     const RPCDetId Rsid = (*detUnitIt).first;
0086     const RPCRoll *roll = dynamic_cast<const RPCRoll *>(rpcGeom->roll(Rsid));
0087
0088     const RPCDigiCollection::Range &range = (*detUnitIt).second;
0089     std::vector<double> sims;
0090     if (allsims.find(Rsid) != allsims.end()) {
0091       sims = allsims[Rsid];
0092     }
0093
0094     int ndigi = 0;
0095     for (RPCDigiCollection::const_iterator digiIt = range.first; digiIt != range.second; ++digiIt) {
0096       StripProf->Fill(digiIt->strip());
0097       BxDist->Fill(digiIt->bx());
0098       // bx for 4 endcaps
0099       if (Rsid.region() == (+1)) {
0100         if (Rsid.station() == 4)
0101           BxDisc_4Plus->Fill(digiIt->bx());
0102       } else if (Rsid.region() == (-1)) {
0103         if (Rsid.station() == 4)
0104           BxDisc_4Min->Fill(digiIt->bx());
0105       }
0106
0107       // Fill timing information
0108       const double digiTime = digiIt->hasTime() ? digiIt->time() : digiIt->bx() * 25;
0109       hDigiTimeAll->Fill(digiTime);
0110       if (digiIt->hasTime()) {
0111         hDigiTime->Fill(digiTime);
0112         if (roll->isIRPC())
0113           hDigiTimeIRPC->Fill(digiTime);
0114         else
0115           hDigiTimeNoIRPC->Fill(digiTime);
0116       }
0117     }
0118
0119     if (sims.size() == 1 && ndigi == 1) {
0120       double dis = roll->centreOfStrip(range.first->strip()).x() - sims[0];
0121       Res->Fill(dis);
0122
0123       if (Rsid.region() == 0) {
0124         if (Rsid.ring() == -2)
0125           ResWmin2->Fill(dis);
0126         else if (Rsid.ring() == -1)
0127           ResWmin1->Fill(dis);
0128         else if (Rsid.ring() == 0)
0129           ResWzer0->Fill(dis);
0130         else if (Rsid.ring() == 1)
0131           ResWplu1->Fill(dis);
0132         else if (Rsid.ring() == 2)
0133           ResWplu2->Fill(dis);
0134         // barrel layers
0135         if (Rsid.station() == 1 && Rsid.layer() == 1)
0136           ResLayer1_barrel->Fill(dis);
0137         if (Rsid.station() == 1 && Rsid.layer() == 2)
0138           ResLayer2_barrel->Fill(dis);
0139         if (Rsid.station() == 2 && Rsid.layer() == 1)
0140           ResLayer3_barrel->Fill(dis);
0141         if (Rsid.station() == 2 && Rsid.layer() == 2)
0142           ResLayer4_barrel->Fill(dis);
0143         if (Rsid.station() == 3)
0144           ResLayer5_barrel->Fill(dis);
0145         if (Rsid.station() == 4)
0146           ResLayer6_barrel->Fill(dis);
0147       }
0148       // endcap layers residuals
0149       if (Rsid.region() != 0) {
0150         if (Rsid.ring() == 2) {
0151           if (abs(Rsid.station()) == 1) {
0152             if (Rsid.roll() == 1)
0153               Res_Endcap1_Ring2_A->Fill(dis);
0154             if (Rsid.roll() == 2)
0155               Res_Endcap1_Ring2_B->Fill(dis);
0156             if (Rsid.roll() == 3)
0157               Res_Endcap1_Ring2_C->Fill(dis);
0158           }
0159           if (abs(Rsid.station()) == 2 || abs(Rsid.station()) == 3) {
0160             if (Rsid.roll() == 1)
0161               Res_Endcap23_Ring2_A->Fill(dis);
0162             if (Rsid.roll() == 2)
0163               Res_Endcap23_Ring2_B->Fill(dis);
0164             if (Rsid.roll() == 3)
0165               Res_Endcap23_Ring2_C->Fill(dis);
0166           }
0167         }
0168         if (Rsid.ring() == 3) {
0169           if (Rsid.roll() == 1)
0170             Res_Endcap123_Ring3_A->Fill(dis);
0171           if (Rsid.roll() == 2)
0172             Res_Endcap123_Ring3_B->Fill(dis);
0173           if (Rsid.roll() == 3)
0174             Res_Endcap123_Ring3_C->Fill(dis);
0175         }
0176       }
0177
0178       if (Rsid.region() == (+1)) {
0179         if (Rsid.station() == 1)
0180           ResDplu1->Fill(dis);
0181         else if (Rsid.station() == 2)
0182           ResDplu2->Fill(dis);
0183         else if (Rsid.station() == 3)
0184           ResDplu3->Fill(dis);
0185         else if (Rsid.station() == 4)
0186           ResDplu4->Fill(dis);
0187       }
0188       if (Rsid.region() == (-1)) {
0189         if (Rsid.station() == 1)
0190           ResDmin1->Fill(dis);
0191         else if (Rsid.station() == 2)
0192           ResDmin2->Fill(dis);
0193         else if (Rsid.station() == 3)
0194           ResDmin3->Fill(dis);
0195         else if (Rsid.station() == 4)
0196           ResDmin4->Fill(dis);
0197       }
0198     }
0199   }
0200 }
0201
0202 void RPCDigiValid::bookHistograms(DQMStore::IBooker &booker, edm::Run const &run, edm::EventSetup const &eSetup) {
0203   booker.setCurrentFolder("RPCDigisV/RPCDigis");
0204
0205   xyview = booker.book2D("X_Vs_Y_View", "X_Vs_Y_View", 155, -775., 775., 155, -775., 775.);
0206
0207   xyvDplu4 = booker.book2D("Dplu4_XvsY", "Dplu4_XvsY", 155, -775., 775., 155, -775., 775.);
0208   xyvDmin4 = booker.book2D("Dmin4_XvsY", "Dmin4_XvsY", 155, -775., 775., 155, -775., 775.);
0209
0210   rzview = booker.book2D("R_Vs_Z_View", "R_Vs_Z_View", 216, -1080., 1080., 52, 260., 780.);
0211   Res = booker.book1D("Digi_SimHit_difference", "Digi_SimHit_difference", 300, -8, 8);
0212   ResWmin2 = booker.book1D("W_Min2_Residuals", "W_Min2_Residuals", 400, -8, 8);
0213   ResWmin1 = booker.book1D("W_Min1_Residuals", "W_Min1_Residuals", 400, -8, 8);
0214   ResWzer0 = booker.book1D("W_Zer0_Residuals", "W_Zer0_Residuals", 400, -8, 8);
0215   ResWplu1 = booker.book1D("W_Plu1_Residuals", "W_Plu1_Residuals", 400, -8, 8);
0216   ResWplu2 = booker.book1D("W_Plu2_Residuals", "W_Plu2_Residuals", 400, -8, 8);
0217
0218   ResLayer1_barrel = booker.book1D("ResLayer1_barrel", "ResLayer1_barrel", 400, -8, 8);
0219   ResLayer2_barrel = booker.book1D("ResLayer2_barrel", "ResLayer2_barrel", 400, -8, 8);
0220   ResLayer3_barrel = booker.book1D("ResLayer3_barrel", "ResLayer3_barrel", 400, -8, 8);
0221   ResLayer4_barrel = booker.book1D("ResLayer4_barrel", "ResLayer4_barrel", 400, -8, 8);
0222   ResLayer5_barrel = booker.book1D("ResLayer5_barrel", "ResLayer5_barrel", 400, -8, 8);
0223   ResLayer6_barrel = booker.book1D("ResLayer6_barrel", "ResLayer6_barrel", 400, -8, 8);
0224
0225   BxDist = booker.book1D("Bunch_Crossing", "Bunch_Crossing", 20, -10., 10.);
0226   StripProf = booker.book1D("Strip_Profile", "Strip_Profile", 100, 0, 100);
0227
0228   BxDisc_4Plus = booker.book1D("BxDisc_4Plus", "BxDisc_4Plus", 20, -10., 10.);
0229   BxDisc_4Min = booker.book1D("BxDisc_4Min", "BxDisc_4Min", 20, -10., 10.);
0230
0231   // endcap residuals
0232   ResDmin1 = booker.book1D("Disk_Min1_Residuals", "Disk_Min1_Residuals", 400, -8, 8);
0233   ResDmin2 = booker.book1D("Disk_Min2_Residuals", "Disk_Min2_Residuals", 400, -8, 8);
0234   ResDmin3 = booker.book1D("Disk_Min3_Residuals", "Disk_Min3_Residuals", 400, -8, 8);
0235   ResDplu1 = booker.book1D("Disk_Plu1_Residuals", "Disk_Plu1_Residuals", 400, -8, 8);
0236   ResDplu2 = booker.book1D("Disk_Plu2_Residuals", "Disk_Plu2_Residuals", 400, -8, 8);
0237   ResDplu3 = booker.book1D("Disk_Plu3_Residuals", "Disk_Plu3_Residuals", 400, -8, 8);
0238
0239   ResDmin4 = booker.book1D("Disk_Min4_Residuals", "Disk_Min4_Residuals", 400, -8, 8);
0240   ResDplu4 = booker.book1D("Disk_Plu4_Residuals", "Disk_Plu4_Residuals", 400, -8, 8);
0241
0242   Res_Endcap1_Ring2_A = booker.book1D("Res_Endcap1_Ring2_A", "Res_Endcap1_Ring2_A", 400, -8, 8);
0243   Res_Endcap1_Ring2_B = booker.book1D("Res_Endcap1_Ring2_B", "Res_Endcap1_Ring2_B", 400, -8, 8);
0244   Res_Endcap1_Ring2_C = booker.book1D("Res_Endcap1_Ring2_C", "Res_Endcap1_Ring2_C", 400, -8, 8);
0245
0246   Res_Endcap23_Ring2_A = booker.book1D("Res_Endcap23_Ring2_A", "Res_Endcap23_Ring2_A", 400, -8, 8);
0247   Res_Endcap23_Ring2_B = booker.book1D("Res_Endcap23_Ring2_B", "Res_Endcap23_Ring2_B", 400, -8, 8);
0248   Res_Endcap23_Ring2_C = booker.book1D("Res_Endcap23_Ring2_C", "Res_Endcap23_Ring2_C", 400, -8, 8);
0249
0250   Res_Endcap123_Ring3_A = booker.book1D("Res_Endcap123_Ring3_A", "Res_Endcap123_Ring3_A", 400, -8, 8);
0251   Res_Endcap123_Ring3_B = booker.book1D("Res_Endcap123_Ring3_B", "Res_Endcap123_Ring3_B", 400, -8, 8);
0252   Res_Endcap123_Ring3_C = booker.book1D("Res_Endcap123_Ring3_C", "Res_Endcap123_Ring3_C", 400, -8, 8);
0253
0254   // Timing informations
0255   hDigiTimeAll =
0256       booker.book1D("DigiTimeAll", "Digi time including present electronics;Digi time (ns)", 100, -12.5, 12.5);
0257   hDigiTime = booker.book1D("DigiTime", "Digi time only with timing information;Digi time (ns)", 100, -12.5, 12.5);
0258   hDigiTimeIRPC = booker.book1D("DigiTimeIRPC", "IRPC Digi time;Digi time (ns)", 100, -12.5, 12.5);
0259   hDigiTimeNoIRPC = booker.book1D("DigiTimeNoIRPC", "non-IRPC Digi time;Digi time (ns)", 100, -12.5, 12.5);
0260 }