File indexing completed on 2023-10-25 10:06:52
0001 #include "Validation/MuonGEMHits/plugins/GEMSimHitValidation.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0004
0005 GEMSimHitValidation::GEMSimHitValidation(const edm::ParameterSet& pset)
0006 : GEMBaseValidation(pset, "GEMSimHitValidation") {
0007 const auto& simhit_pset = pset.getParameterSet("gemSimHit");
0008 const auto& simhit_tag = simhit_pset.getParameter<edm::InputTag>("inputTag");
0009 simhit_token_ = consumes<edm::PSimHitContainer>(simhit_tag);
0010
0011 tof_range_ = pset.getUntrackedParameter<std::vector<Double_t> >("TOFRange");
0012 geomToken_ = esConsumes<GEMGeometry, MuonGeometryRecord>();
0013 geomTokenBeginRun_ = esConsumes<GEMGeometry, MuonGeometryRecord, edm::Transition::BeginRun>();
0014 }
0015
0016 GEMSimHitValidation::~GEMSimHitValidation() {}
0017
0018 void GEMSimHitValidation::bookHistograms(DQMStore::IBooker& booker, edm::Run const& run, edm::EventSetup const& setup) {
0019 const GEMGeometry* gem = &setup.getData(geomTokenBeginRun_);
0020
0021
0022 booker.setCurrentFolder("GEM/SimHits");
0023
0024 TString tof_xtitle = "Time of flight [ns]";
0025 TString tof_ytitle = "Entries";
0026
0027 for (const auto& region : gem->regions()) {
0028 Int_t region_id = region->region();
0029
0030 for (const auto& station : region->stations()) {
0031 Int_t station_id = station->station();
0032
0033 const auto [tof_min, tof_max] = getTOFRange(station_id);
0034 ME2IdsKey key2{region_id, station_id};
0035
0036 me_tof_mu_[key2] = bookHist1D(
0037 booker, key2, "sim_tof_muon", "SimHit TOF (Muon only)", 20, tof_min, tof_max, tof_xtitle, tof_ytitle);
0038
0039 me_tof_others_[key2] = bookHist1D(
0040 booker, key2, "sim_tof_others", "SimHit TOF (Other Particles)", 20, tof_min, tof_max, tof_xtitle, tof_ytitle);
0041 }
0042 }
0043
0044 if (detail_plot_) {
0045 for (const auto& region : gem->regions()) {
0046 Int_t region_id = region->region();
0047
0048 for (const auto& station : region->stations()) {
0049 Int_t station_id = station->station();
0050
0051 const auto [tof_min, tof_max] = getTOFRange(station_id);
0052 const auto& superChamberVec = station->superChambers();
0053 if (superChamberVec.empty() || superChamberVec.front() == nullptr) {
0054 edm::LogError(kLogCategory_) << "Super chambers missing or null for region = " << region_id
0055 << " and station = " << station_id;
0056 } else {
0057 const GEMSuperChamber* super_chamber = superChamberVec.front();
0058
0059 for (const auto& chamber : super_chamber->chambers()) {
0060 Int_t layer_id = chamber->id().layer();
0061 ME3IdsKey key3{region_id, station_id, layer_id};
0062
0063 me_detail_tof_[key3] = bookHist1D(
0064 booker, key3, "sim_tof", "Time of Flight of Muon SimHits", 40, tof_min, tof_max, tof_xtitle, tof_ytitle);
0065
0066 me_detail_tof_mu_[key3] = bookHist1D(
0067 booker, key3, "sim_tof_muon", "SimHit TOF (Muon only)", 40, tof_min, tof_max, tof_xtitle, tof_ytitle);
0068 }
0069 }
0070 }
0071 }
0072 }
0073
0074
0075 TString eloss_xtitle = "Energy loss [eV]";
0076 TString eloss_ytitle = "Entries / 0.5 keV";
0077
0078
0079
0080 for (const auto& region : gem->regions()) {
0081 for (const auto& station : region->stations()) {
0082 Int_t station_id = station->station();
0083 if (me_eloss_mu_.find(station_id) != me_eloss_mu_.end())
0084 continue;
0085
0086 auto eloss_mu_name = TString::Format("sim_eloss_muon_GE%d1", station_id);
0087 auto eloss_mu_title = TString::Format("SimHit Energy Loss (Muon only) : GE%d1", station_id);
0088
0089 me_eloss_mu_[station_id] =
0090 booker.book1D(eloss_mu_name, eloss_mu_title + ";" + eloss_xtitle + ";" + eloss_ytitle, 20, 0.0, 10.0);
0091
0092 auto eloss_others_name = TString::Format("sim_eloss_others_GE%d1", station_id);
0093 auto eloss_others_title = TString::Format("SimHit Energy Loss (Other Particles) : GE%d1", station_id);
0094
0095 me_eloss_others_[station_id] =
0096 booker.book1D(eloss_others_name, eloss_others_title + ";" + eloss_xtitle + ";" + eloss_ytitle, 20, 0.0, 10.0);
0097 }
0098 }
0099
0100 if (detail_plot_) {
0101 for (const auto& region : gem->regions()) {
0102 Int_t region_id = region->region();
0103 for (const auto& station : region->stations()) {
0104 Int_t station_id = station->station();
0105 const auto& superChamberVec = station->superChambers();
0106 if (superChamberVec.empty() || superChamberVec.front() == nullptr) {
0107 edm::LogError(kLogCategory_) << "Super chambers missing or null for region = " << region_id
0108 << " and station = " << station_id;
0109 } else {
0110 for (const auto& chamber : superChamberVec.front()->chambers()) {
0111 Int_t layer_id = chamber->id().layer();
0112 ME3IdsKey key3{region_id, station_id, layer_id};
0113
0114 me_detail_eloss_[key3] = bookHist1D(
0115 booker, key3, "sim_eloss", "SimHit Energy Loss", 60, 0.0, 6000.0, eloss_xtitle, eloss_ytitle);
0116
0117 me_detail_eloss_mu_[key3] = bookHist1D(booker,
0118 key3,
0119 "sim_eloss_muon",
0120 "SimHit Energy Loss (Muon Only)",
0121 60,
0122 0.0,
0123 6000.0,
0124 eloss_xtitle,
0125 eloss_ytitle);
0126
0127 }
0128 }
0129 }
0130 }
0131 }
0132
0133
0134 for (const auto& region : gem->regions()) {
0135 Int_t region_id = region->region();
0136
0137 if (detail_plot_)
0138 me_detail_occ_zr_[region_id] = bookZROccupancy(booker, region_id, "sim_simhit", "SimHit");
0139
0140 for (const auto& station : region->stations()) {
0141 Int_t station_id = station->station();
0142 ME2IdsKey key2{region_id, station_id};
0143
0144 if (detail_plot_) {
0145 me_detail_occ_det_[key2] = bookDetectorOccupancy(booker, key2, station, "sim_simhit", "SimHit");
0146 me_detail_occ_det_mu_[key2] = bookDetectorOccupancy(booker, key2, station, "sim_muon_simhit", "Muon SimHit");
0147 }
0148
0149 const auto& superChamberVec = station->superChambers();
0150 if (superChamberVec.empty() || superChamberVec.front() == nullptr) {
0151 edm::LogError(kLogCategory_) << "Super chambers missing or null for region = " << region_id
0152 << " and station = " << station_id;
0153 } else {
0154 const GEMSuperChamber* super_chamber = superChamberVec.front();
0155 for (const auto& chamber : super_chamber->chambers()) {
0156 Int_t layer_id = chamber->id().layer();
0157 ME3IdsKey key3{region_id, station_id, layer_id};
0158
0159 me_occ_eta_mu_[key3] = bookHist1D(booker,
0160 key3,
0161 "sim_muon_occ_eta",
0162 "Muon SimHit Eta Occupancy",
0163 16,
0164 eta_range_[station_id * 2 + 0],
0165 eta_range_[station_id * 2 + 1],
0166 "#eta");
0167
0168 me_occ_phi_mu_[key3] =
0169 bookHist1D(booker, key3, "sim_muon_occ_phi", "Muon SimHit Phi Occupancy", 36, -5, 355, "#phi [degrees]");
0170
0171 me_occ_pid_[key3] = bookPIDHist(booker, key3, "sim_occ_pid", "Particle population");
0172
0173 if (detail_plot_)
0174 me_detail_occ_xy_[key3] = bookXYOccupancy(booker, key3, "sim_simhit", "SimHit");
0175 }
0176 }
0177 }
0178 }
0179 }
0180
0181 std::tuple<Double_t, Double_t> GEMSimHitValidation::getTOFRange(Int_t station_id) {
0182 UInt_t start_index = station_id * 2;
0183 Double_t tof_min = tof_range_[start_index];
0184 Double_t tof_max = tof_range_[start_index + 1];
0185 return std::make_tuple(tof_min, tof_max);
0186 }
0187
0188 void GEMSimHitValidation::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0189 const GEMGeometry* gem = &setup.getData(geomToken_);
0190
0191 edm::Handle<edm::PSimHitContainer> simhit_container;
0192 event.getByToken(simhit_token_, simhit_container);
0193 if (not simhit_container.isValid()) {
0194 edm::LogError(kLogCategory_) << "Cannot get GEMHits by Token simhitLabel" << std::endl;
0195 return;
0196 }
0197
0198 for (const auto& simhit : *simhit_container.product()) {
0199 const GEMDetId gemid(simhit.detUnitId());
0200
0201 if (gem->idToDet(gemid) == nullptr) {
0202 edm::LogError(kLogCategory_) << "SimHit did not matched with GEM Geometry." << std::endl;
0203 continue;
0204 }
0205
0206 Int_t region_id = gemid.region();
0207 Int_t station_id = gemid.station();
0208 Int_t layer_id = gemid.layer();
0209 Int_t chamber_id = gemid.chamber();
0210 Int_t ieta = gemid.ieta();
0211 Int_t num_layers = gemid.nlayers();
0212
0213 ME2IdsKey key2{region_id, station_id};
0214 ME3IdsKey key3{region_id, station_id, layer_id};
0215
0216 GlobalPoint&& simhit_global_pos = gem->idToDet(gemid)->surface().toGlobal(simhit.localPosition());
0217
0218 Float_t simhit_g_r = simhit_global_pos.perp();
0219 Float_t simhit_g_x = simhit_global_pos.x();
0220 Float_t simhit_g_y = simhit_global_pos.y();
0221 Float_t simhit_g_abs_z = std::fabs(simhit_global_pos.z());
0222 Float_t simhit_g_eta = std::fabs(simhit_global_pos.eta());
0223 Float_t simhit_g_phi = toDegree(simhit_global_pos.phi());
0224
0225 Float_t energy_loss = kEnergyCF_ * simhit.energyLoss();
0226 energy_loss = energy_loss > 10 ? 9.9 : energy_loss;
0227 Float_t tof = simhit.timeOfFlight();
0228 Int_t pid = simhit.particleType();
0229 Int_t pid_idx = getPidIdx(pid);
0230
0231
0232 Int_t bin_x = getDetOccBinX(num_layers, chamber_id, layer_id);
0233
0234 bool is_muon_simhit = isMuonSimHit(simhit);
0235 if (is_muon_simhit) {
0236 me_tof_mu_[key2]->Fill(tof);
0237 me_eloss_mu_[station_id]->Fill(energy_loss);
0238 me_occ_eta_mu_[key3]->Fill(simhit_g_eta);
0239 me_occ_phi_mu_[key3]->Fill(simhit_g_phi);
0240 } else {
0241 me_tof_others_[key2]->Fill(tof);
0242 me_eloss_others_[station_id]->Fill(energy_loss);
0243 }
0244
0245 me_occ_pid_[key3]->Fill(pid_idx);
0246
0247 if (detail_plot_) {
0248 me_detail_tof_[key3]->Fill(tof);
0249 me_detail_eloss_[key3]->Fill(energy_loss);
0250
0251 me_detail_occ_zr_[region_id]->Fill(simhit_g_abs_z, simhit_g_r);
0252 me_detail_occ_det_[key2]->Fill(bin_x, ieta);
0253 me_detail_occ_xy_[key3]->Fill(simhit_g_x, simhit_g_y);
0254
0255 if (is_muon_simhit) {
0256 me_detail_tof_mu_[key3]->Fill(tof);
0257 me_detail_eloss_mu_[key3]->Fill(energy_loss);
0258 me_detail_occ_det_mu_[key2]->Fill(bin_x, ieta);
0259 }
0260
0261 }
0262 }
0263 }
