Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-27 18:27:36

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   // NOTE Occupancy
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   // NOTE Occupancy
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       }  // end loop over layer ids
0115     }    // end loop over station ids
0116   }      // end loop over region ids
0117 
0118   // NOTE Bunch Crossing
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         }  // chamber loop
0143       }    // station loop
0144     }      // region loop
0145   }        // detail plot
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       // bunch crossing
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       }  // detail_plot_
0247     }
0248   }  // end loop over range iters
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   // NOTE
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   }  // simhit_container loop
0331 }