Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-02 23:43:48

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   // NOTE Time of flight
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     }  // station loop
0042   }    // region loop
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           }  // chamber loop
0069         }    // end else
0070       }      // station loop
0071     }        // region loop
0072   }          // detail plot
0073 
0074   // NOTE Energy Loss
0075   TString eloss_xtitle = "Energy loss [eV]";
0076   TString eloss_ytitle = "Entries / 0.5 keV";
0077 
0078   // Demonstrator chamber means we could have a missing station,
0079   // process both regions to make sure we include it
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     }  // station loop
0098   }    // region loop
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           }  // chamber loop
0128         }    // end else
0129       }      // station loop
0130     }        // region loop
0131   }          // detail plot
0132 
0133   // NOTE Occupancy
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         }  // layer loop
0176       }    // end else
0177     }      // station loop
0178   }        // region loop
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     // NOTE Fill MonitorElement
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     }  // detail_plot
0262   }    // simhit loop
0263 }