Back to home page

Project CMSSW displayed by LXR

 
 

    


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   // NOTE Bunch Crossing
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         }  // chamber loop
0069       }    // station loop
0070     }      // region loop
0071   }
0072 
0073   // NOTE Occupancy
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         }  // chamber
0154       }    // end else
0155     }      // station looop
0156   }        // region loop
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   // NOTE
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   }  // range loop
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   // NOTE
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   }  // simhit_container loop
0313 }