File indexing completed on 2023-03-17 11:28:08
0001 #include "Validation/MuonGEMDigis/plugins/GEMPadDigiClusterValidation.h"
0002 #include <TMath.h>
0003
0004 GEMPadDigiClusterValidation::GEMPadDigiClusterValidation(const edm::ParameterSet& pset)
0005 : GEMBaseValidation(pset, "GEMPadDigiClusterValidation") {
0006 const auto& pad_cluster_pset = pset.getParameterSet("gemPadCluster");
0007 const auto& pad_cluster_tag = pad_cluster_pset.getParameter<edm::InputTag>("inputTag");
0008
0009 const auto& simhit_pset = pset.getParameterSet("gemSimHit");
0010 const auto& simhit_tag = simhit_pset.getParameter<edm::InputTag>("inputTag");
0011 simhit_token_ = consumes<edm::PSimHitContainer>(simhit_tag);
0012
0013 const auto& digisimlink_tag = pset.getParameter<edm::InputTag>("gemDigiSimLink");
0014 digisimlink_token_ = consumes<edm::DetSetVector<GEMDigiSimLink>>(digisimlink_tag);
0015
0016 pad_cluster_token_ = consumes<GEMPadDigiClusterCollection>(pad_cluster_tag);
0017 geomToken_ = esConsumes<GEMGeometry, MuonGeometryRecord>();
0018 geomTokenBeginRun_ = esConsumes<GEMGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0019 }
0020
0021 void GEMPadDigiClusterValidation::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 booker.setCurrentFolder("GEM/PadCluster");
0032
0033 TString cls_title = "Cluster Size Distribution";
0034 TString cls_x_title = "Cluster size";
0035
0036 me_cls_ = booker.book1D("cls", cls_title + ";" + cls_x_title + ";" + "Entries", 10, 0, 10);
0037
0038
0039 for (const auto& region : gem->regions()) {
0040 Int_t region_id = region->region();
0041
0042 if (detail_plot_)
0043 me_detail_occ_zr_.emplace(region_id, bookZROccupancy(booker, region_id, "pad", "Pad Cluster"));
0044
0045 for (const auto& station : region->stations()) {
0046 Int_t station_id = station->station();
0047 ME2IdsKey key2{region_id, station_id};
0048
0049 if (detail_plot_) {
0050 me_detail_occ_det_[key2] = bookDetectorOccupancy(booker, key2, station, "pad", "Pad Cluster");
0051 me_detail_pad_cluster_occ_det_[key2] =
0052 bookDetectorOccupancy(booker, key2, station, "sim_matched", "Pad Cluster");
0053 }
0054
0055 const auto& superChamberVec = station->superChambers();
0056 if (superChamberVec.empty()) {
0057 edm::LogError(kLogCategory_) << "Super chambers missing for region = " << region_id
0058 << " and station = " << station_id;
0059 continue;
0060 }
0061 const GEMSuperChamber* super_chamber = superChamberVec.front();
0062 if (super_chamber == nullptr) {
0063 edm::LogError(kLogCategory_) << "Failed to find super chamber for region = " << region_id
0064 << " and station = " << station_id;
0065 continue;
0066 }
0067 for (const auto& chamber : super_chamber->chambers()) {
0068 Int_t layer_id = chamber->id().layer();
0069 ME3IdsKey key3{region_id, station_id, layer_id};
0070
0071 const auto& etaPartitionVec = chamber->etaPartitions();
0072 if (etaPartitionVec.empty() || etaPartitionVec.front() == nullptr) {
0073 edm::LogError(kLogCategory_)
0074 << "Eta partition missing or null for region, station, super chamber, chamber = (" << region_id << ", "
0075 << station_id << ", " << super_chamber->id() << ", " << chamber->id() << ")";
0076 continue;
0077 }
0078 Int_t num_pads = etaPartitionVec.front()->npads();
0079
0080 me_total_cluster_[key3] =
0081 bookHist1D(booker, key3, "total_pad_cluster", "Number of pad digi cluster per event", 20, 0, 20);
0082
0083 me_pad_cluster_occ_eta_[key3] = bookHist1D(booker,
0084 key3,
0085 "sim_matched_occ_eta",
0086 "Matched Pad Cluster Eta Occupancy",
0087 16,
0088 eta_range_[station_id * 2 + 0],
0089 eta_range_[station_id * 2 + 1],
0090 "#eta");
0091
0092 me_pad_cluster_occ_phi_[key3] = bookHist1D(
0093 booker, key3, "sim_matched_occ_phi", "Matched Pad Cluster Phi Occupancy", 36, -5, 355, "#phi [degrees]");
0094
0095 if (detail_plot_) {
0096 me_detail_occ_xy_[key3] = bookXYOccupancy(booker, key3, "pad", "Pad Cluster");
0097
0098 me_detail_occ_phi_pad_[key3] = bookHist2D(booker,
0099 key3,
0100 "occ_phi_pad",
0101 "Pad Cluster Occupancy",
0102 280,
0103 -M_PI,
0104 M_PI,
0105 num_pads / 2,
0106 0,
0107 num_pads,
0108 "#phi [rad]",
0109 "Pad number");
0110
0111 me_detail_occ_pad_[key3] =
0112 bookHist1D(booker, key3, "occ_pad", "Pad Cluster Occupancy", num_pads, 0, num_pads, "Pad number");
0113 }
0114 }
0115 }
0116 }
0117
0118
0119 if (detail_plot_) {
0120 for (const auto& region : gem->regions()) {
0121 Int_t region_id = region->region();
0122 for (const auto& station : region->stations()) {
0123 Int_t station_id = station->station();
0124
0125 const auto& superChamberVec = station->superChambers();
0126 if (superChamberVec.empty()) {
0127 edm::LogError(kLogCategory_) << "Super chambers missing for region = " << region_id
0128 << " and station = " << station_id;
0129 continue;
0130 }
0131 const GEMSuperChamber* super_chamber = superChamberVec.front();
0132 if (super_chamber == nullptr) {
0133 edm::LogError(kLogCategory_) << "Failed to find super chamber for region = " << region_id
0134 << " and station = " << station_id;
0135 continue;
0136 }
0137 for (const auto& chamber : super_chamber->chambers()) {
0138 Int_t layer_id = chamber->id().layer();
0139 ME3IdsKey key3(region_id, station_id, layer_id);
0140 me_detail_bx_[key3] =
0141 bookHist1D(booker, key3, "bx", "Pad Cluster Bunch Crossing", 5, -2, 3, "Bunch crossing");
0142 }
0143 }
0144 }
0145 }
0146 }
0147
0148 GEMPadDigiClusterValidation::~GEMPadDigiClusterValidation() {}
0149
0150 Bool_t GEMPadDigiClusterValidation::matchClusterAgainstSimHit(GEMPadDigiClusterCollection::const_iterator cluster,
0151 Int_t simhit_pad) {
0152 for (auto pad : cluster->pads()) {
0153 if (pad == simhit_pad) {
0154 return true;
0155 }
0156 }
0157 return false;
0158 }
0159
0160 void GEMPadDigiClusterValidation::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0161 const auto& gemH = setup.getHandle(geomToken_);
0162 if (!gemH.isValid()) {
0163 edm::LogError(kLogCategory_) << "Failed to initialize GEM geometry.";
0164 return;
0165 }
0166 const GEMGeometry* gem = gemH.product();
0167
0168 edm::Handle<GEMPadDigiClusterCollection> collection;
0169 event.getByToken(pad_cluster_token_, collection);
0170 if (not collection.isValid()) {
0171 edm::LogError(kLogCategory_) << "Cannot get pads by label GEMPadToken.";
0172 return;
0173 }
0174
0175 edm::Handle<edm::DetSetVector<GEMDigiSimLink>> digiSimLink;
0176 event.getByToken(digisimlink_token_, digiSimLink);
0177 if (not digiSimLink.isValid()) {
0178 edm::LogError(kLogCategory_) << "Failed to get GEMDigiSimLink." << std::endl;
0179 return;
0180 }
0181
0182 edm::Handle<edm::PSimHitContainer> simhit_container;
0183 event.getByToken(simhit_token_, simhit_container);
0184 if (not simhit_container.isValid()) {
0185 edm::LogError(kLogCategory_) << "Failed to get PSimHitContainer." << std::endl;
0186 return;
0187 }
0188
0189 std::map<ME3IdsKey, Int_t> total_cluster;
0190 for (auto range_iter = collection->begin(); range_iter != collection->end(); range_iter++) {
0191 GEMDetId gemid = (*range_iter).first;
0192 const auto& range = (*range_iter).second;
0193
0194 if (gem->idToDet(gemid) == nullptr) {
0195 edm::LogError(kLogCategory_) << "Getting DetId failed. Discard this gem pad hit. "
0196 << "Maybe it comes from unmatched geometry." << std::endl;
0197 continue;
0198 }
0199
0200 const GEMEtaPartition* roll = gem->etaPartition(gemid);
0201 const BoundPlane& surface = roll->surface();
0202
0203 Int_t region_id = gemid.region();
0204 Int_t station_id = gemid.station();
0205 Int_t layer_id = gemid.layer();
0206 Int_t chamber_id = gemid.chamber();
0207 Int_t ieta = gemid.ieta();
0208 Int_t num_layers = gemid.nlayers();
0209
0210 ME2IdsKey key2(region_id, station_id);
0211 ME3IdsKey key3(region_id, station_id, layer_id);
0212
0213 for (auto digi = range.first; digi != range.second; ++digi) {
0214 const auto& padsVec = digi->pads();
0215 if (padsVec.empty()) {
0216 edm::LogError(kLogCategory_) << "Pads missing for digi from GEM ID = " << gemid;
0217 continue;
0218 }
0219 Int_t pad = padsVec[0];
0220
0221 total_cluster[key3]++;
0222
0223
0224 Int_t bx = digi->bx();
0225 Int_t cls = digi->pads().size();
0226
0227 const LocalPoint& local_pos = roll->centreOfPad(pad);
0228 const GlobalPoint& global_pos = surface.toGlobal(local_pos);
0229
0230 Float_t g_r = global_pos.perp();
0231 Float_t g_phi = global_pos.phi();
0232 Float_t g_x = global_pos.x();
0233 Float_t g_y = global_pos.y();
0234 Float_t g_abs_z = std::fabs(global_pos.z());
0235
0236 Int_t bin_x = getDetOccBinX(num_layers, chamber_id, layer_id);
0237
0238 me_cls_->Fill(cls);
0239 if (detail_plot_) {
0240 me_detail_occ_zr_[region_id]->Fill(g_abs_z, g_r);
0241 me_detail_occ_det_[key2]->Fill(bin_x, ieta);
0242 me_detail_occ_xy_[key3]->Fill(g_x, g_y);
0243 me_detail_occ_phi_pad_[key3]->Fill(g_phi, pad);
0244 me_detail_occ_pad_[key3]->Fill(pad);
0245 me_detail_bx_[key3]->Fill(bx);
0246 }
0247 }
0248 }
0249
0250 for (const auto& region : gem->regions()) {
0251 Int_t region_id = region->region();
0252 for (const auto& station : region->stations()) {
0253 Int_t station_id = station->station();
0254 const auto& superChamberVec = station->superChambers();
0255 if (superChamberVec.empty()) {
0256 edm::LogError(kLogCategory_) << "Super chambers missing for region = " << region_id
0257 << " and station = " << station_id;
0258 continue;
0259 }
0260 const GEMSuperChamber* super_chamber = superChamberVec.front();
0261 if (super_chamber == nullptr) {
0262 edm::LogError(kLogCategory_) << "Failed to find super chamber for region = " << region_id
0263 << " and station = " << station_id;
0264 continue;
0265 }
0266 for (const auto& chamber : super_chamber->chambers()) {
0267 Int_t layer_id = chamber->id().layer();
0268 ME3IdsKey key3{region_id, station_id, layer_id};
0269 me_total_cluster_[key3]->Fill(total_cluster[key3]);
0270 }
0271 }
0272 }
0273
0274
0275 for (const auto& simhit : *simhit_container.product()) {
0276 if (not isMuonSimHit(simhit))
0277 continue;
0278 if (gem->idToDet(simhit.detUnitId()) == nullptr) {
0279 edm::LogError(kLogCategory_) << "SimHit did not match with GEMGeometry." << std::endl;
0280 continue;
0281 }
0282
0283 GEMDetId simhit_gemid(simhit.detUnitId());
0284
0285 Int_t region_id = simhit_gemid.region();
0286 Int_t station_id = simhit_gemid.station();
0287 Int_t layer_id = simhit_gemid.layer();
0288 Int_t chamber_id = simhit_gemid.chamber();
0289 Int_t ieta = simhit_gemid.ieta();
0290 Int_t num_layers = simhit_gemid.nlayers();
0291
0292 ME2IdsKey key2{region_id, station_id};
0293 ME3IdsKey key3{region_id, station_id, layer_id};
0294
0295 const GEMEtaPartition* roll = gem->etaPartition(simhit_gemid);
0296
0297 const auto& simhit_local_pos = simhit.localPosition();
0298 const auto& simhit_global_pos = roll->surface().toGlobal(simhit_local_pos);
0299
0300 Float_t simhit_g_eta = std::abs(simhit_global_pos.eta());
0301 Float_t simhit_g_phi = toDegree(simhit_global_pos.phi());
0302
0303 auto simhit_trackId = simhit.trackId();
0304
0305 Int_t bin_x = getDetOccBinX(num_layers, chamber_id, layer_id);
0306
0307 auto links = digiSimLink->find(simhit_gemid);
0308 if (links == digiSimLink->end())
0309 continue;
0310
0311 Int_t simhit_strip = -1;
0312 for (const auto& link : *links) {
0313 if (simhit_trackId == link.getTrackId()) {
0314 simhit_strip = link.getStrip();
0315 break;
0316 }
0317 }
0318 Int_t simhit_pad = roll->padOfStrip(simhit_strip);
0319 auto range = collection->get(simhit_gemid);
0320 for (auto cluster = range.first; cluster != range.second; ++cluster) {
0321 if (matchClusterAgainstSimHit(cluster, simhit_pad)) {
0322 me_pad_cluster_occ_eta_[key3]->Fill(simhit_g_eta);
0323 me_pad_cluster_occ_phi_[key3]->Fill(simhit_g_phi);
0324 if (detail_plot_) {
0325 me_detail_pad_cluster_occ_det_[key2]->Fill(bin_x, ieta);
0326 }
0327 break;
0328 }
0329 }
0330 }
0331 }