File indexing completed on 2023-10-25 10:06:52
0001 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0002 #include "Validation/MuonGEMDigis/plugins/GEMStripDigiValidation.h"
0003
0004 GEMStripDigiValidation::GEMStripDigiValidation(const edm::ParameterSet& pset)
0005 : GEMBaseValidation(pset, "GEMStripDigiValidation") {
0006 const auto& strip_pset = pset.getParameterSet("gemStripDigi");
0007 const auto& strip_tag = strip_pset.getParameter<edm::InputTag>("inputTag");
0008 strip_token_ = consumes<GEMDigiCollection>(strip_tag);
0009
0010 const auto& simhit_pset = pset.getParameterSet("gemSimHit");
0011 const auto& simhit_tag = simhit_pset.getParameter<edm::InputTag>("inputTag");
0012 simhit_token_ = consumes<edm::PSimHitContainer>(simhit_tag);
0013
0014 const auto& digisimlink_tag = pset.getParameter<edm::InputTag>("gemDigiSimLink");
0015 digisimlink_token_ = consumes<edm::DetSetVector<GEMDigiSimLink>>(digisimlink_tag);
0016
0017 geomToken_ = esConsumes<GEMGeometry, MuonGeometryRecord>();
0018 geomTokenBeginRun_ = esConsumes<GEMGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0019 }
0020
0021 void GEMStripDigiValidation::bookHistograms(DQMStore::IBooker& booker,
0022 edm::Run const& run,
0023 edm::EventSetup const& setup) {
0024 const auto gemH = setup.getHandle(geomTokenBeginRun_);
0025 if (!gemH.isValid()) {
0026 edm::LogError(kLogCategory_) << "Failed to initialize GEM geometry.";
0027 return;
0028 }
0029 const GEMGeometry* gem = gemH.product();
0030
0031
0032 booker.setCurrentFolder("GEM/Digis");
0033
0034 if (detail_plot_) {
0035 me_detail_bx_ = booker.book1D("bx", "Strip Digi Bunch Crossing", 5, -2.5, 2.5);
0036
0037 for (const auto& region : gem->regions()) {
0038 if (region == nullptr) {
0039 edm::LogError(kLogCategory_) << "Null region";
0040 continue;
0041 }
0042 Int_t region_id = region->region();
0043 for (const auto& station : region->stations()) {
0044 if (station == nullptr) {
0045 edm::LogError(kLogCategory_) << "Null station for region = " << region_id;
0046 continue;
0047 }
0048 Int_t station_id = station->station();
0049
0050 const auto& superChamberVec = station->superChambers();
0051 if (superChamberVec.empty()) {
0052 edm::LogError(kLogCategory_) << "Super chambers missing for region = " << region_id
0053 << " and station = " << station_id;
0054 continue;
0055 }
0056 const GEMSuperChamber* super_chamber = superChamberVec.front();
0057 if (super_chamber == nullptr) {
0058 edm::LogError(kLogCategory_) << "Failed to find super chamber for region = " << region_id
0059 << " and station = " << station_id;
0060 continue;
0061 }
0062 for (const auto& chamber : super_chamber->chambers()) {
0063 Int_t layer_id = chamber->id().layer();
0064 ME3IdsKey key3(region_id, station_id, layer_id);
0065
0066 me_detail_bx_layer_[key3] =
0067 bookHist1D(booker, key3, "bx", "Strip Digi Bunch Crossing", 5, -2.5, 2.5, "Bunch crossing");
0068 }
0069 }
0070 }
0071 }
0072
0073
0074 if (detail_plot_)
0075 me_detail_total_strip_all_ =
0076 booker.book1D("total_strips_per_event", "Number of strip digi per event", 50, -0.5, 395.5);
0077
0078 for (const auto& region : gem->regions()) {
0079 Int_t region_id = region->region();
0080
0081 if (detail_plot_)
0082 me_detail_occ_zr_[region_id] = bookZROccupancy(booker, region_id, "strip", "Strip Digi");
0083
0084 for (const auto& station : region->stations()) {
0085 Int_t station_id = station->station();
0086 ME2IdsKey key2{region_id, station_id};
0087
0088 me_occ_pid_[key2] = bookPIDHist(booker, key2, "sim_occ_pid", "Particle population");
0089
0090 if (detail_plot_) {
0091 me_detail_total_strip_[key2] =
0092 bookHist1D(booker, key2, "total_strips_per_event", "Number of strip digs per event", 50, -0.5, 99.5);
0093
0094 me_detail_occ_det_[key2] = bookDetectorOccupancy(booker, key2, station, "strip", "Strip Digi");
0095
0096 me_detail_strip_occ_det_[key2] =
0097 bookDetectorOccupancy(booker, key2, station, "sim_matched_strip", "Matched Strip Digi");
0098 }
0099
0100 const auto& superChamberVec = station->superChambers();
0101 if (superChamberVec.empty() || superChamberVec[0] == nullptr) {
0102 edm::LogError(kLogCategory_) << "Super chambers missing or null for region = " << region_id
0103 << " and station = " << station_id;
0104 } else {
0105 for (const auto& etaPart : superChamberVec[0]->chambers()[0]->etaPartitions()) {
0106 Int_t ieta = etaPart->id().ieta();
0107 ME3IdsKey key{region_id, station_id, ieta};
0108 me_occ_pid_eta_[key] = bookPIDHist(booker, key2, ieta, "sim_occ_pid", "Particle population");
0109 }
0110 for (const auto& chamber : superChamberVec[0]->chambers()) {
0111 if (chamber == nullptr) {
0112 edm::LogError(kLogCategory_) << "Null chamber for region, station, super chamber = (" << region_id << ", "
0113 << station_id << ", " << superChamberVec[0]->id() << ")";
0114 continue;
0115 }
0116 Int_t layer_id = chamber->id().layer();
0117 ME3IdsKey key3{region_id, station_id, layer_id};
0118
0119 const auto& etaPartitionsVec = chamber->etaPartitions();
0120 if (etaPartitionsVec.empty() || etaPartitionsVec.front() == nullptr) {
0121 edm::LogError(kLogCategory_)
0122 << "Eta partition missing or null for region, station, super chamber, chamber = (" << region_id << ", "
0123 << station_id << ", " << superChamberVec[0]->id() << ", " << chamber->id() << ")";
0124 continue;
0125 }
0126
0127 if (detail_plot_) {
0128 Int_t num_strips = etaPartitionsVec.front()->nstrips();
0129
0130 me_detail_strip_occ_eta_[key3] = bookHist1D(booker,
0131 key3,
0132 "sim_matched_occ_eta",
0133 "Matched Strip Eta Occupancy",
0134 16,
0135 eta_range_[station_id * 2 + 0],
0136 eta_range_[station_id * 2 + 1],
0137 "#eta");
0138
0139 me_detail_strip_occ_phi_[key3] = bookHist1D(
0140 booker, key3, "sim_matched_occ_phi", "Matched Strip Phi Occupancy", 36, -5, 355, "#phi [degrees]");
0141
0142 me_detail_occ_xy_[key3] = bookXYOccupancy(booker, key3, "strip", "Strip Digi");
0143
0144 me_detail_occ_strip_[key3] = bookHist1D(booker,
0145 key3,
0146 "occ_strip",
0147 "Strip Digi Occupancy per strip number",
0148 num_strips,
0149 0.5,
0150 num_strips + 0.5,
0151 "strip number");
0152 }
0153 }
0154 }
0155 }
0156 }
0157 }
0158
0159 GEMStripDigiValidation::~GEMStripDigiValidation() {}
0160
0161 void GEMStripDigiValidation::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0162 const auto& gemH = setup.getHandle(geomToken_);
0163 if (!gemH.isValid()) {
0164 edm::LogError(kLogCategory_) << "Failed to initialize GEM geometry.";
0165 return;
0166 }
0167 const GEMGeometry* gem = gemH.product();
0168
0169 edm::Handle<edm::DetSetVector<GEMDigiSimLink>> digiSimLink;
0170 event.getByToken(digisimlink_token_, digiSimLink);
0171 if (not digiSimLink.isValid()) {
0172 edm::LogError(kLogCategory_) << "Failed to get GEMDigiSimLink." << std::endl;
0173 return;
0174 }
0175
0176 edm::Handle<edm::PSimHitContainer> simhit_container;
0177 event.getByToken(simhit_token_, simhit_container);
0178 if (not simhit_container.isValid()) {
0179 edm::LogError(kLogCategory_) << "Failed to get PSimHitContainer." << std::endl;
0180 return;
0181 }
0182
0183 edm::Handle<GEMDigiCollection> digi_collection;
0184 event.getByToken(strip_token_, digi_collection);
0185 if (not digi_collection.isValid()) {
0186 edm::LogError(kLogCategory_) << "Cannot get strips by Token stripToken." << std::endl;
0187 return;
0188 }
0189
0190
0191 Int_t total_strip = 0;
0192 std::map<ME2IdsKey, Int_t> total_strip_2IdsMap;
0193 for (const auto digi_pair : *digi_collection) {
0194 GEMDetId id = digi_pair.first;
0195 if (gem->idToDet(id) == nullptr) {
0196 edm::LogError(kLogCategory_) << "Getting DetId failed. Discard this gem strip hit. Maybe it comes "
0197 << "from unmatched geometry." << std::endl;
0198 continue;
0199 }
0200
0201 Int_t region_id = id.region();
0202 Int_t layer_id = id.layer();
0203 Int_t station_id = id.station();
0204 Int_t chamber_id = id.chamber();
0205 Int_t ieta = id.ieta();
0206 Int_t num_layers = id.nlayers();
0207
0208 ME2IdsKey key2{region_id, station_id};
0209 ME3IdsKey key3{region_id, station_id, layer_id};
0210 Int_t bin_x = getDetOccBinX(num_layers, chamber_id, layer_id);
0211
0212 const BoundPlane& surface = gem->idToDet(id)->surface();
0213 const GEMEtaPartition* roll = gem->etaPartition(id);
0214
0215 const GEMDigiCollection::Range& range = digi_pair.second;
0216 auto links = digiSimLink->find(id);
0217 if (links == digiSimLink->end())
0218 continue;
0219
0220 for (auto digi = range.first; digi != range.second; ++digi) {
0221 total_strip++;
0222 total_strip_2IdsMap[key2]++;
0223 Int_t strip = digi->strip();
0224 Int_t bx = digi->bx();
0225 bx = bx < -10 ? -10 : bx;
0226 bx = bx > 10 ? 10 : bx;
0227
0228 GlobalPoint strip_global_pos = surface.toGlobal(roll->centreOfStrip(strip));
0229
0230 Float_t digi_g_x = strip_global_pos.x();
0231 Float_t digi_g_y = strip_global_pos.y();
0232 Float_t digi_g_r = strip_global_pos.perp();
0233 Float_t digi_g_abs_z = std::abs(strip_global_pos.z());
0234
0235 if (detail_plot_) {
0236 me_detail_bx_->Fill(bx);
0237 me_detail_bx_layer_[key3]->Fill(bx);
0238
0239 me_detail_occ_zr_[region_id]->Fill(digi_g_abs_z, digi_g_r);
0240 me_detail_occ_det_[key2]->Fill(bin_x, ieta);
0241 me_detail_occ_xy_[key3]->Fill(digi_g_x, digi_g_y);
0242 me_detail_occ_strip_[key3]->Fill(strip);
0243 }
0244 }
0245 }
0246 if (detail_plot_) {
0247 me_detail_total_strip_all_->Fill(total_strip);
0248 for (const auto& region : gem->regions()) {
0249 Int_t region_id = region->region();
0250 for (const auto& station : region->stations()) {
0251 Int_t station_id = station->station();
0252 ME2IdsKey key2{region_id, station_id};
0253 me_detail_total_strip_[key2]->Fill(total_strip_2IdsMap[key2]);
0254 }
0255 }
0256 }
0257
0258
0259 for (const auto& simhit : *simhit_container.product()) {
0260 if (gem->idToDet(simhit.detUnitId()) == nullptr) {
0261 edm::LogError(kLogCategory_) << "SimHit did not match with GEMGeometry." << std::endl;
0262 continue;
0263 }
0264
0265 GEMDetId simhit_gemid(simhit.detUnitId());
0266
0267 Int_t region_id = simhit_gemid.region();
0268 Int_t station_id = simhit_gemid.station();
0269 Int_t layer_id = simhit_gemid.layer();
0270 Int_t chamber_id = simhit_gemid.chamber();
0271 Int_t ieta = simhit_gemid.ieta();
0272 Int_t num_layers = simhit_gemid.nlayers();
0273
0274 ME3IdsKey key{region_id, station_id, ieta};
0275 ME2IdsKey key2{region_id, station_id};
0276 ME3IdsKey key3{region_id, station_id, layer_id};
0277
0278 const GEMEtaPartition* roll = gem->etaPartition(simhit_gemid);
0279
0280 const auto& simhit_local_pos = simhit.localPosition();
0281 const auto& simhit_global_pos = roll->surface().toGlobal(simhit_local_pos);
0282
0283 Float_t simhit_g_eta = std::abs(simhit_global_pos.eta());
0284 Float_t simhit_g_phi = toDegree(simhit_global_pos.phi());
0285
0286 auto simhit_trackId = simhit.trackId();
0287
0288 Int_t bin_x = getDetOccBinX(num_layers, chamber_id, layer_id);
0289
0290 auto links = digiSimLink->find(simhit_gemid);
0291 if (links == digiSimLink->end())
0292 continue;
0293
0294 for (const auto& link : *links) {
0295 if (simhit_trackId == link.getTrackId()) {
0296 Int_t pid = simhit.particleType();
0297 Int_t pid_idx = getPidIdx(pid);
0298
0299 me_occ_pid_[key2]->Fill(pid_idx);
0300 me_occ_pid_eta_[key]->Fill(pid_idx);
0301
0302 if (detail_plot_) {
0303 if (isMuonSimHit(simhit)) {
0304 me_detail_strip_occ_eta_[key3]->Fill(simhit_g_eta);
0305 me_detail_strip_occ_phi_[key3]->Fill(simhit_g_phi);
0306 me_detail_strip_occ_det_[key2]->Fill(bin_x, ieta);
0307 }
0308 }
0309 break;
0310 }
0311 }
0312 }
0313 }
