MuonHits/src/ME0SimHitMatcher.cc

0001 #include "Validation/MuonHits/interface/ME0SimHitMatcher.h"
0002
0003 using namespace std;
0004
0005 ME0SimHitMatcher::ME0SimHitMatcher(const edm::ParameterSet& ps, edm::ConsumesCollector&& iC)
0006     : MuonSimHitMatcher(ps, std::move(iC)) {
0007   simHitPSet_ = ps.getParameterSet("me0SimHit");
0008   verbose_ = simHitPSet_.getParameter<int>("verbose");
0009   simMuOnly_ = simHitPSet_.getParameter<bool>("simMuOnly");
0010   discardEleHits_ = simHitPSet_.getParameter<bool>("discardEleHits");
0011
0012   simHitInput_ = iC.consumes<edm::PSimHitContainer>(simHitPSet_.getParameter<edm::InputTag>("inputTag"));
0013   geomToken_ = iC.esConsumes<ME0Geometry, MuonGeometryRecord>();
0014 }
0015
0016 /// initialize the event
0017 void ME0SimHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0018   geometry_ = &iSetup.getData(geomToken_);
0019   MuonSimHitMatcher::init(iEvent, iSetup);
0020 }
0021
0022 /// do the matching
0023 void ME0SimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) {
0024   // instantiates the track ids and simhits
0025   MuonSimHitMatcher::match(track, vertex);
0026
0027   matchSimHitsToSimTrack();
0028
0029   if (verbose_) {
0030     edm::LogInfo("ME0SimHitMatcher") << "nTrackIds " << track_ids_.size() << " nSelectedME0SimHits " << hits_.size();
0031     edm::LogInfo("ME0SimHitMatcher") << "detids ME0 " << detIds().size();
0032
0033     const auto& me0_ch_ids = detIds();
0034     for (const auto& id : me0_ch_ids) {
0035       const auto& me0_simhits = MuonSimHitMatcher::hitsInChamber(id);
0036       const auto& me0_simhits_gp = simHitsMeanPosition(me0_simhits);
0037       edm::LogInfo("ME0SimHitMatcher") << "me0chid " << ME0DetId(id) << ": nHits " << me0_simhits.size() << " phi "
0038                                        << me0_simhits_gp.phi() << " nCh " << chamber_to_hits_[id].size();
0039       const auto& strips = hitStripsInDetId(id);
0040       edm::LogInfo("ME0SimHitMatcher") << "nStrip " << strips.size();
0041       edm::LogInfo("ME0SimHitMatcher") << "strips : ";
0042       for (const auto& p : strips) {
0043         edm::LogInfo("ME0SimHitMatcher") << p;
0044       }
0045     }
0046   }
0047 }
0048
0049 void ME0SimHitMatcher::matchSimHitsToSimTrack() {
0050   for (const auto& track_id : track_ids_) {
0051     for (const auto& h : simHits_) {
0052       if (h.trackId() != track_id)
0053         continue;
0054       int pdgid = h.particleType();
0055       if (simMuOnly_ && std::abs(pdgid) != 13)
0056         continue;
0057       // discard electron hits in the ME0 chambers
0058       if (discardEleHits_ && std::abs(pdgid) == 11)
0059         continue;
0060
0061       const ME0DetId& layer_id(h.detUnitId());
0062       detid_to_hits_[h.detUnitId()].push_back(h);
0063       hits_.push_back(h);
0064       chamber_to_hits_[layer_id.layerId().rawId()].push_back(h);
0065       superChamber_to_hits_[layer_id.chamberId().rawId()].push_back(h);
0066     }
0067   }
0068
0069   // find pads with hits
0070   const auto& detids = detIds();
0071   for (const auto& d : detids) {
0072     ME0DetId id(d);
0073     const auto& hits = hitsInDetId(d);
0074     const auto& roll = dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(id);
0075     // int max_npads = roll->npads();
0076     set<int> pads;
0077     for (const auto& h : hits) {
0078       const LocalPoint& lp = h.entryPoint();
0079       pads.insert(1 + static_cast<int>(roll->padTopology().channel(lp)));
0080     }
0081     detids_to_pads_[d] = pads;
0082   }
0083 }
0084
0085 std::set<unsigned int> ME0SimHitMatcher::detIds() const {
0086   std::set<unsigned int> result;
0087   for (const auto& p : detid_to_hits_)
0088     result.insert(p.first);
0089   return result;
0090 }
0091
0092 std::set<unsigned int> ME0SimHitMatcher::chamberIds() const {
0093   std::set<unsigned int> result;
0094   for (const auto& p : chamber_to_hits_)
0095     result.insert(p.first);
0096   return result;
0097 }
0098
0099 std::set<unsigned int> ME0SimHitMatcher::superChamberIds() const {
0100   std::set<unsigned int> result;
0101   for (const auto& p : superChamber_to_hits_)
0102     result.insert(p.first);
0103   return result;
0104 }
0105
0106 const edm::PSimHitContainer& ME0SimHitMatcher::hitsInSuperChamber(unsigned int detid) const {
0107   if (superChamber_to_hits_.find(detid) == superChamber_to_hits_.end())
0108     return no_hits_;
0109   return superChamber_to_hits_.at(detid);
0110 }
0111
0112 int ME0SimHitMatcher::nLayersWithHitsInSuperChamber(unsigned int detid) const {
0113   set<int> layers_with_hits;
0114   const auto& hits = hitsInSuperChamber(detid);
0115   for (const auto& h : hits) {
0116     const ME0DetId& idd(h.detUnitId());
0117     layers_with_hits.insert(idd.layer());
0118   }
0119   return layers_with_hits.size();
0120 }
0121
0122 std::set<unsigned int> ME0SimHitMatcher::superChamberIdsCoincidences(int min_n_layers) const {
0123   set<unsigned int> result;
0124   // loop over the super chamber Ids
0125   for (const auto& p : superChamberIds()) {
0126     if (nLayersWithHitsInSuperChamber(p) >= min_n_layers) {
0127       result.insert(p);
0128     }
0129   }
0130   return result;
0131 }
0132
0133 int ME0SimHitMatcher::nCoincidenceChambers(int min_n_layers) const {
0134   return superChamberIdsCoincidences(min_n_layers).size();
0135 }
0136
0137 float ME0SimHitMatcher::simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const {
0138   if (sim_hits.empty())
0139     return -1.f;
0140
0141   float sums = 0.f;
0142   size_t n = 0;
0143   for (const auto& h : sim_hits) {
0144     const LocalPoint& lp = h.entryPoint();
0145     float s;
0146     const auto& d = h.detUnitId();
0147     s = dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(d)->strip(lp);
0148     sums += s;
0149     ++n;
0150   }
0151   if (n == 0)
0152     return -1.f;
0153   return sums / n;
0154 }
0155
0156 std::set<int> ME0SimHitMatcher::hitStripsInDetId(unsigned int detid, int margin_n_strips) const {
0157   set<int> result;
0158   const auto& simhits = hitsInDetId(detid);
0159   ME0DetId id(detid);
0160   int max_nstrips = dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(id)->nstrips();
0161   for (const auto& h : simhits) {
0162     const LocalPoint& lp = h.entryPoint();
0163     int central_strip =
0164         1 + static_cast<int>(dynamic_cast<const ME0Geometry*>(geometry_)->etaPartition(id)->topology().channel(lp));
0165     int smin = central_strip - margin_n_strips;
0166     smin = (smin > 0) ? smin : 1;
0167     int smax = central_strip + margin_n_strips;
0168     smax = (smax <= max_nstrips) ? smax : max_nstrips;
0169     for (int ss = smin; ss <= smax; ++ss)
0170       result.insert(ss);
0171   }
0172   return result;
0173 }
0174
0175 std::set<int> ME0SimHitMatcher::hitPadsInDetId(unsigned int detid) const {
0176   set<int> none;
0177   if (detids_to_pads_.find(detid) == detids_to_pads_.end())
0178     return none;
0179   return detids_to_pads_.at(detid);
0180 }
0181
0182 std::set<int> ME0SimHitMatcher::hitPartitions() const {
0183   std::set<int> result;
0184
0185   const auto& detids = detIds();
0186   for (const auto& id : detids) {
0187     ME0DetId idd(id);
0188     result.insert(idd.roll());
0189   }
0190   return result;
0191 }
0192
0193 int ME0SimHitMatcher::nPadsWithHits() const {
0194   int result = 0;
0195   const auto& pad_ids = detIds();
0196   for (const auto& id : pad_ids) {
0197     result += hitPadsInDetId(id).size();
0198   }
0199   return result;
0200 }