MuonGEMRecHits/plugins/GEMRecHitValidation.cc

0001 #include "Validation/MuonGEMRecHits/plugins/GEMRecHitValidation.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0005
0006 GEMRecHitValidation::GEMRecHitValidation(const edm::ParameterSet& pset)
0007     : GEMBaseValidation(pset, "GEMRecHitsValidation") {
0008   const auto& rechit_pset = pset.getParameterSet("gemRecHit");
0009   const auto& rechit_tag = rechit_pset.getParameter<edm::InputTag>("inputTag");
0010   rechit_token_ = consumes<GEMRecHitCollection>(rechit_tag);
0011
0012   const auto& simhit_pset = pset.getParameterSet("gemSimHit");
0013   const auto& simhit_tag = simhit_pset.getParameter<edm::InputTag>("inputTag");
0014
0015   const auto& digisimlink_tag = pset.getParameter<edm::InputTag>("gemDigiSimLink");
0016   digisimlink_token_ = consumes<edm::DetSetVector<GEMDigiSimLink>>(digisimlink_tag);
0017
0018   simhit_token_ = consumes<edm::PSimHitContainer>(simhit_tag);
0019   geomToken_ = esConsumes<GEMGeometry, MuonGeometryRecord>();
0020   geomTokenBeginRun_ = esConsumes<GEMGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0021 }
0022
0023 GEMRecHitValidation::~GEMRecHitValidation() {}
0024
0025 void GEMRecHitValidation::bookHistograms(DQMStore::IBooker& booker, edm::Run const& run, edm::EventSetup const& setup) {
0026   const GEMGeometry* gem = &setup.getData(geomTokenBeginRun_);
0027
0028   // NOTE Cluster Size
0029   booker.setCurrentFolder("GEM/RecHits");
0030
0031   TString cls_title = "Cluster Size Distribution";
0032   TString cls_x_title = "Cluster size";
0033
0034   if (detail_plot_) {
0035     me_detail_cls_total_ = booker.book1D("cls", cls_title + ";" + cls_x_title + ";" + "Entries", 10, 0.5, 10.5);
0036
0037     for (const auto& station : gem->regions()[0]->stations()) {
0038       Int_t station_id = station->station();
0039       for (const auto& roll : station->superChambers()[0]->chambers()[0]->etaPartitions()) {
0040         Int_t ieta = roll->id().ieta();
0041         ME2IdsKey key{station_id, ieta};
0042         me_detail_cls_roll_[key] = booker.book1D(Form("cls_GE%d1-E%d", station_id, ieta),
0043                                                  Form("Cluster Size Distribution : GE%d1-E%d", station_id, ieta),
0044                                                  10,
0045                                                  0.5,
0046                                                  10.5);
0047       }
0048     }
0049
0050     for (const auto& region : gem->regions()) {
0051       Int_t region_id = region->region();
0052
0053       for (const auto& station : region->stations()) {
0054         Int_t station_id = station->station();
0055
0056         const auto& superChamberVec = station->superChambers();
0057         if (superChamberVec.empty() || superChamberVec[0] == nullptr) {
0058           edm::LogError(kLogCategory_) << "Super chambers missing or null for region = " << region_id
0059                                        << " and station = " << station_id;
0060         } else {
0061           for (const auto& chamber : superChamberVec[0]->chambers()) {
0062             Int_t layer_id = chamber->id().layer();
0063
0064             for (const auto& roll : chamber->etaPartitions()) {
0065               Int_t ieta = roll->id().ieta();
0066               ME4IdsKey key4{region_id, station_id, layer_id, ieta};
0067
0068               me_detail_cls_[key4] = bookHist1D(booker, key4, "cls", "Cluster Size Distribution", 11, -0.5, 10.5);
0069             }  // roll loop
0070           }    // chamber loop
0071         }      // end else
0072       }        // station loop
0073     }          // region loop
0074   }            // detail plot
0075
0076   // NOTE Residual
0077   for (const auto& station : gem->regions()[0]->stations()) {
0078     Int_t station_id = station->station();
0079     for (const auto& roll : station->superChambers()[0]->chambers()[0]->etaPartitions()) {
0080       Int_t ieta = roll->id().ieta();
0081       ME2IdsKey key{station_id, ieta};
0082
0083       me_residual_y_[key] = booker.book1D(Form("residual_y_GE%d1-E%d", station_id, ieta),
0084                                           Form("Residual in Y : GE%d1-E%d; Residual in Y [cm]", station_id, ieta),
0085                                           60,
0086                                           -15,
0087                                           15);
0088
0089       me_residual_rphi_[key] =
0090           booker.book1D(Form("residual_rphi_GE%d1-E%d", station_id, ieta),
0091                         Form("Residual in R #times #phi : GE%d1-E%d; Residual in r #times #phi [cm]", station_id, ieta),
0092                         60,
0093                         -15,
0094                         15);
0095     }
0096   }
0097
0098   if (detail_plot_) {
0099     for (const auto& region : gem->regions()) {
0100       Int_t region_id = region->region();
0101
0102       for (const auto& station : region->stations()) {
0103         Int_t station_id = station->station();
0104
0105         const auto& superChamberVec = station->superChambers();
0106         if (!superChamberVec.empty() && superChamberVec[0] != nullptr) {
0107           for (const auto& chamber : superChamberVec[0]->chambers()) {
0108             Int_t layer_id = chamber->id().layer();
0109
0110             for (const auto& roll : chamber->etaPartitions()) {
0111               Int_t ieta = roll->id().ieta();
0112               ME4IdsKey key4{region_id, station_id, layer_id, ieta};
0113
0114               me_detail_residual_y_[key4] =
0115                   bookHist1D(booker, key4, "residual_y", "Residual in y", 60, -15, 15, "Residual in y [cm]");
0116
0117               me_detail_residual_rphi_[key4] = bookHist1D(booker,
0118                                                           key4,
0119                                                           "residual_rphi",
0120                                                           "Residual in r #times #phi",
0121                                                           60,
0122                                                           -15,
0123                                                           15,
0124                                                           "Residual in r #times #phi [cm]");
0125             }  // roll loop
0126           }    // chamber loop
0127         }      // end if
0128       }        // station loop
0129     }          // region loop
0130   }            // detail plot
0131
0132   // NOTE Pull
0133   if (detail_plot_) {
0134     for (const auto& station : gem->regions()[0]->stations()) {
0135       Int_t station_id = station->station();
0136       for (const auto& roll : station->superChambers()[0]->chambers()[0]->etaPartitions()) {
0137         Int_t ieta = roll->id().ieta();
0138         ME2IdsKey key{station_id, ieta};
0139
0140         me_detail_pull_x_[key] = booker.book1D(
0141             Form("pull_x_GE%d1-E%d", station_id, ieta), Form("Pull in X : GE%d1-E%d", station_id, ieta), 60, -3, 3);
0142
0143         me_detail_pull_y_[key] = booker.book1D(
0144             Form("pull_y_GE%d1-E%d", station_id, ieta), Form("Pull in Y : GE%d1-E%d", station_id, ieta), 60, -3, 3);
0145       }
0146     }
0147
0148     for (const auto& region : gem->regions()) {
0149       Int_t region_id = region->region();
0150
0151       for (const auto& station : region->stations()) {
0152         Int_t station_id = station->station();
0153
0154         const auto& superChamberVec = station->superChambers();
0155         if (!superChamberVec.empty() && superChamberVec[0] != nullptr) {
0156           for (const auto& chamber : superChamberVec[0]->chambers()) {
0157             Int_t layer_id = chamber->id().layer();
0158
0159             for (const auto& roll : chamber->etaPartitions()) {
0160               Int_t ieta = roll->id().ieta();
0161               ME4IdsKey key4{region_id, station_id, layer_id, ieta};
0162
0163               me_detail_pull_x_la_[key4] = bookHist1D(booker, key4, "pull_x", "Pull in x", 60, -3, 3);
0164
0165               me_detail_pull_y_la_[key4] = bookHist1D(booker, key4, "pull_y", "Pull in y", 60, -3, 3);
0166             }  // roll loop
0167           }    // chamber loop
0168         }      // end if
0169       }        // station loop
0170     }          // region loop
0171   }            // detail plot
0172
0173   // NOTE Occupancy
0174   for (const auto& region : gem->regions()) {
0175     Int_t region_id = region->region();
0176
0177     if (detail_plot_)
0178       me_detail_occ_zr_[region_id] = bookZROccupancy(booker, region_id, "rechit", "RecHit");
0179
0180     for (const auto& station : region->stations()) {
0181       Int_t station_id = station->station();
0182       ME2IdsKey key2{region_id, station_id};
0183
0184       if (detail_plot_)
0185         me_detail_rechit_occ_det_[key2] = bookDetectorOccupancy(booker, key2, station, "sim_matched", "Matched RecHit");
0186
0187       const auto& superChamberVec = station->superChambers();
0188       if (!superChamberVec.empty() && superChamberVec[0] != nullptr) {
0189         for (const auto& chamber : superChamberVec[0]->chambers()) {
0190           Int_t layer_id = chamber->id().layer();
0191           ME3IdsKey key3{region_id, station_id, layer_id};
0192
0193           Int_t num_eta_partitions = chamber->nEtaPartitions();
0194
0195           me_rechit_occ_eta_[key3] = bookHist1D(booker,
0196                                                 key3,
0197                                                 "sim_matched_occ_eta",
0198                                                 "Matched RecHit Eta Occupancy",
0199                                                 16,
0200                                                 eta_range_[station_id * 2 + 0],
0201                                                 eta_range_[station_id * 2 + 1],
0202                                                 "|#eta|");
0203
0204           me_rechit_occ_phi_[key3] =
0205               bookHist1D(booker, key3, "sim_matched_occ_phi", "Matched RecHit Phi Occupancy", 36, -5, 355, "#phi");
0206
0207           if (detail_plot_) {
0208             me_detail_total_rechit_[key3] =
0209                 bookHist1D(booker, key3, "total_rechit", "Number of rec hits per event", 25, -0.5, 24.5);
0210
0211             me_detail_occ_pid_[key3] = bookPIDHist(booker, key3, "sim_occ_pid", "Number of entreis for each particle");
0212
0213             me_detail_occ_ieta_[key3] = bookHist1D(booker,
0214                                                    key3,
0215                                                    "occ_ieta",
0216                                                    "Rechit Occupancy per eta partition",
0217                                                    num_eta_partitions,
0218                                                    0.5,
0219                                                    num_eta_partitions + 0.5);
0220
0221             me_detail_occ_phi_[key3] = bookHist1D(booker, key3, "occ_phi", "Rechit Phi Occupancy", 108, -5, 355);
0222
0223             me_detail_occ_xy_[key3] = bookXYOccupancy(booker, key3, "rechit", "RecHit");
0224
0225             me_detail_occ_polar_[key3] = bookPolarOccupancy(booker, key3, "rechit", "RecHit");
0226           }
0227         }  // chamber loop
0228       }    // end if
0229     }      // station loop
0230   }        // region_loop
0231 }
0232
0233 Bool_t GEMRecHitValidation::matchRecHitAgainstSimHit(GEMRecHitCollection::const_iterator rechit, Int_t simhit_strip) {
0234   Int_t cls = rechit->clusterSize();
0235   Int_t rechit_first_strip = rechit->firstClusterStrip();
0236
0237   if (cls == 1) {
0238     return simhit_strip == rechit_first_strip;
0239   } else {
0240     Int_t rechit_last_strip = rechit_first_strip + cls - 1;
0241     return (simhit_strip >= rechit_first_strip) and (simhit_strip <= rechit_last_strip);
0242   }
0243 }
0244
0245 void GEMRecHitValidation::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0246   const GEMGeometry* gem = &setup.getData(geomToken_);
0247
0248   edm::Handle<edm::DetSetVector<GEMDigiSimLink>> digiSimLink;
0249   event.getByToken(digisimlink_token_, digiSimLink);
0250   if (not digiSimLink.isValid()) {
0251     edm::LogError(kLogCategory_) << "Failed to get GEMDigiSimLink." << std::endl;
0252     return;
0253   }
0254
0255   edm::Handle<edm::PSimHitContainer> simhit_container;
0256   event.getByToken(simhit_token_, simhit_container);
0257   if (not simhit_container.isValid()) {
0258     edm::LogError(kLogCategory_) << "Failed to get PSimHitContainer." << std::endl;
0259     return;
0260   }
0261
0262   edm::Handle<GEMRecHitCollection> rechit_collection;
0263   event.getByToken(rechit_token_, rechit_collection);
0264   if (not rechit_collection.isValid()) {
0265     edm::LogError(kLogCategory_) << "Failed to get GEMRecHitCollection" << std::endl;
0266     return;
0267   }
0268
0269   std::map<ME3IdsKey, Int_t> total_rechit;
0270   for (const auto& rechit : *rechit_collection) {
0271     GEMDetId gem_id{rechit.gemId()};
0272     Int_t region_id = gem_id.region();
0273     Int_t station_id = gem_id.station();
0274     Int_t layer_id = gem_id.layer();
0275     Int_t ieta = gem_id.ieta();
0276     ME4IdsKey key4{region_id, station_id, layer_id, ieta};
0277
0278     ME2IdsKey key{station_id, ieta};
0279     ME2IdsKey key2{region_id, station_id};
0280     ME3IdsKey key3{region_id, station_id, layer_id};
0281
0282     const BoundPlane& surface = gem->idToDet(gem_id)->surface();
0283     GlobalPoint&& rechit_global_pos = surface.toGlobal(rechit.localPosition());
0284
0285     Float_t rechit_g_x = rechit_global_pos.x();
0286     Float_t rechit_g_y = rechit_global_pos.y();
0287     Float_t rechit_g_abs_z = std::fabs(rechit_global_pos.z());
0288     Float_t rechit_g_r = rechit_global_pos.perp();
0289     Float_t rechit_g_phi = toDegree(rechit_global_pos.phi());
0290
0291     Int_t first_strip = rechit.firstClusterStrip();
0292     Int_t cls = rechit.clusterSize();
0293     cls = cls > 10 ? 10 : cls;
0294
0295     total_rechit[key3]++;
0296
0297     if (detail_plot_) {
0298       me_detail_cls_total_->Fill(cls);
0299       me_detail_cls_[key4]->Fill(cls);
0300       me_detail_cls_roll_[key]->Fill(cls);
0301       me_detail_occ_ieta_[key3]->Fill(ieta);
0302       me_detail_occ_phi_[key3]->Fill(rechit_g_phi);
0303       me_detail_occ_zr_[region_id]->Fill(rechit_g_abs_z, rechit_g_r);
0304       me_detail_occ_xy_[key3]->Fill(rechit_g_x, rechit_g_y);
0305       me_detail_occ_polar_[key3]->Fill(rechit_g_phi, rechit_g_r);
0306     }  // detail plot
0307
0308     auto links = digiSimLink->find(gem_id);
0309
0310     if (links == digiSimLink->end())
0311       continue;
0312     std::map<Int_t, Int_t> pid_count;
0313
0314     for (Int_t strip = first_strip; strip < first_strip + cls; strip++) {
0315       for (const auto& link : *links) {
0316         Int_t link_strip = link.getStrip();
0317         if (link_strip == strip) {
0318           Int_t pid = link.getParticleType();
0319           pid_count[pid]++;
0320           break;
0321         }
0322       }
0323     }
0324     if (detail_plot_) {
0325       Int_t max_pid = 0;
0326       Int_t max_count = 0;
0327       for (auto& [pid, count] : pid_count) {
0328         if (max_count < count) {
0329           max_pid = pid;
0330           max_count = count;
0331         }
0332       }
0333       Int_t pid_idx = getPidIdx(max_pid);
0334       me_detail_occ_pid_[key3]->Fill(pid_idx);
0335     }
0336   }
0337
0338   if (detail_plot_) {
0339     for (auto [key, num_total_rechit] : total_rechit) {
0340       me_detail_total_rechit_[key]->Fill(num_total_rechit);
0341     }
0342   }
0343
0344   // NOTE
0345   for (const auto& simhit : *simhit_container.product()) {
0346     if (gem->idToDet(simhit.detUnitId()) == nullptr) {
0347       edm::LogError(kLogCategory_) << "MuonGEMHit didn't matched with GEMGeometry." << std::endl;
0348       continue;
0349     }
0350
0351     GEMDetId simhit_gemid{simhit.detUnitId()};
0352     const BoundPlane& surface = gem->idToDet(simhit_gemid)->surface();
0353
0354     Int_t region_id = simhit_gemid.region();
0355     Int_t station_id = simhit_gemid.station();
0356     Int_t layer_id = simhit_gemid.layer();
0357     Int_t chamber_id = simhit_gemid.chamber();
0358     Int_t ieta = simhit_gemid.ieta();
0359     Int_t num_layers = simhit_gemid.nlayers();
0360
0361     ME2IdsKey key{station_id, ieta};
0362     ME2IdsKey key2{region_id, station_id};
0363     ME3IdsKey key3{region_id, station_id, layer_id};
0364     ME4IdsKey key4{region_id, station_id, layer_id, ieta};
0365
0366     const LocalPoint& simhit_local_pos = simhit.localPosition();
0367     const GlobalPoint& simhit_global_pos = surface.toGlobal(simhit_local_pos);
0368
0369     Float_t simhit_g_abs_eta = std::fabs(simhit_global_pos.eta());
0370     Float_t simhit_g_phi = toDegree(simhit_global_pos.phi());
0371
0372     Int_t det_occ_bin_x = getDetOccBinX(num_layers, chamber_id, layer_id);
0373
0374     auto simhit_trackId = simhit.trackId();
0375
0376     auto links = digiSimLink->find(simhit_gemid);
0377     if (links == digiSimLink->end())
0378       continue;
0379
0380     Int_t simhit_strip = -1;
0381     for (const auto& link : *links) {
0382       if (simhit_trackId == link.getTrackId()) {
0383         simhit_strip = link.getStrip();
0384         break;
0385       }
0386     }
0387
0388     GEMRecHitCollection::range range = rechit_collection->get(simhit_gemid);
0389     for (auto rechit = range.first; rechit != range.second; ++rechit) {
0390       if (gem->idToDet(rechit->gemId()) == nullptr) {
0391         edm::LogError(kLogCategory_) << "GEMRecHit didn't matched with GEMGeometry." << std::endl;
0392         continue;
0393       }
0394
0395       if (not isMuonSimHit(simhit))
0396         continue;
0397
0398       if (matchRecHitAgainstSimHit(rechit, simhit_strip)) {
0399         const LocalPoint& rechit_local_pos = rechit->localPosition();
0400
0401         Float_t resolution_x = std::sqrt(rechit->localPositionError().xx());
0402         Float_t resolution_y = std::sqrt(rechit->localPositionError().yy());
0403
0404         Float_t residual_x = rechit_local_pos.x() - simhit_local_pos.x();
0405         Float_t residual_y = rechit_local_pos.y() - simhit_local_pos.y();
0406         Float_t residual_r = sqrt(pow(residual_x, 2) + pow(residual_y, 2));
0407         Float_t residual_phi = rechit_local_pos.phi() - simhit_local_pos.phi();
0408         Float_t residual_rphi = residual_r * residual_phi;
0409
0410         Float_t pull_x = residual_x / resolution_x;
0411         Float_t pull_y = residual_y / resolution_y;
0412
0413         me_residual_y_[key]->Fill(residual_y);
0414         me_residual_rphi_[key]->Fill(residual_rphi);
0415
0416         me_rechit_occ_eta_[key3]->Fill(simhit_g_abs_eta);
0417         me_rechit_occ_phi_[key3]->Fill(simhit_g_phi);
0418
0419         if (detail_plot_) {
0420           me_detail_rechit_occ_det_[key2]->Fill(det_occ_bin_x, ieta);
0421
0422           me_detail_residual_y_[key4]->Fill(residual_y);
0423           me_detail_residual_rphi_[key4]->Fill(residual_rphi);
0424
0425           me_detail_pull_x_[key]->Fill(pull_x);
0426           me_detail_pull_y_[key]->Fill(pull_y);
0427           me_detail_pull_x_la_[key4]->Fill(pull_x);
0428           me_detail_pull_y_la_[key4]->Fill(pull_y);
0429         }  // detail_plot
0430
0431         // If we find GEMRecHit that matches PSimHit, then exit
0432         // GEMRecHitCollection loop.
0433         break;
0434
0435       }  // if rechit matches against simhit
0436     }    // rechit loop
0437   }      // simhit loop
0438 }