MuonHits/src/DTSimHitMatcher.cc

0001 #include "Validation/MuonHits/interface/DTSimHitMatcher.h"
0002
0003 using namespace std;
0004
0005 DTSimHitMatcher::DTSimHitMatcher(const edm::ParameterSet& ps, edm::ConsumesCollector&& iC)
0006     : MuonSimHitMatcher(ps, std::move(iC)) {
0007   simHitPSet_ = ps.getParameterSet("dtSimHit");
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<DTGeometry, MuonGeometryRecord>();
0014 }
0015
0016 /// initialize the event
0017 void DTSimHitMatcher::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 DTSimHitMatcher::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("DTSimHitMatcher") << "nTrackIds " << track_ids_.size() << " nSelectedDTSimHits " << hits_.size();
0031     edm::LogInfo("DTSimHitMatcher") << "detids DT " << detIds(0).size();
0032
0033     const auto& dt_det_ids = detIds(0);
0034     for (const auto& id : dt_det_ids) {
0035       const auto& dt_simhits = MuonSimHitMatcher::hitsInDetId(id);
0036       const auto& dt_simhits_gp = simHitsMeanPosition(dt_simhits);
0037       edm::LogInfo("DTSimHitMatcher") << "DTWireId " << DTWireId(id) << ": nHits " << dt_simhits.size() << " eta "
0038                                       << dt_simhits_gp.eta() << " phi " << dt_simhits_gp.phi() << " nCh "
0039                                       << chamber_to_hits_[id].size();
0040     }
0041   }
0042 }
0043
0044 void DTSimHitMatcher::matchSimHitsToSimTrack() {
0045   for (const auto& track_id : track_ids_) {
0046     for (const auto& h : simHits_) {
0047       if (h.trackId() != track_id)
0048         continue;
0049       int pdgid = h.particleType();
0050       if (simMuOnly_ && std::abs(pdgid) != 13)
0051         continue;
0052       // discard electron hits in the DT chambers
0053       if (discardEleHits_ && std::abs(pdgid) == 11)
0054         continue;
0055
0056       const DTWireId layer_id(h.detUnitId());
0057       detid_to_hits_[h.detUnitId()].push_back(h);
0058       hits_.push_back(h);
0059       layer_to_hits_[layer_id.layerId().rawId()].push_back(h);
0060       superlayer_to_hits_[layer_id.superlayerId().rawId()].push_back(h);
0061       chamber_to_hits_[layer_id.chamberId().rawId()].push_back(h);
0062     }
0063   }
0064 }
0065
0066 std::set<unsigned int> DTSimHitMatcher::detIds(int dt_type) const {
0067   std::set<unsigned int> result;
0068   for (const auto& p : detid_to_hits_) {
0069     const auto& id = p.first;
0070     if (dt_type > 0) {
0071       DTWireId detId(id);
0072       if (MuonHitHelper::toDTType(detId.wheel(), detId.station()) != dt_type)
0073         continue;
0074     }
0075     result.insert(id);
0076   }
0077   return result;
0078 }
0079
0080 std::set<unsigned int> DTSimHitMatcher::chamberIds(int dt_type) const {
0081   std::set<unsigned int> result;
0082   for (const auto& p : chamber_to_hits_) {
0083     const auto& id = p.first;
0084     if (dt_type > 0) {
0085       DTChamberId detId(id);
0086       if (MuonHitHelper::toDTType(detId.wheel(), detId.station()) != dt_type)
0087         continue;
0088     }
0089     result.insert(id);
0090   }
0091   return result;
0092 }
0093
0094 std::set<unsigned int> DTSimHitMatcher::layerIds() const {
0095   std::set<unsigned int> result;
0096   for (const auto& p : layer_to_hits_)
0097     result.insert(p.first);
0098   return result;
0099 }
0100
0101 std::set<unsigned int> DTSimHitMatcher::superlayerIds() const {
0102   std::set<unsigned int> result;
0103   for (const auto& p : superlayer_to_hits_)
0104     result.insert(p.first);
0105   return result;
0106 }
0107
0108 const edm::PSimHitContainer& DTSimHitMatcher::hitsInLayer(unsigned int detid) const {
0109   if (!MuonHitHelper::isDT(detid))
0110     return no_hits_;
0111
0112   const DTWireId id(detid);
0113   if (layer_to_hits_.find(id.layerId().rawId()) == layer_to_hits_.end())
0114     return no_hits_;
0115   return layer_to_hits_.at(id.layerId().rawId());
0116 }
0117
0118 const edm::PSimHitContainer& DTSimHitMatcher::hitsInSuperLayer(unsigned int detid) const {
0119   if (!MuonHitHelper::isDT(detid))
0120     return no_hits_;
0121
0122   const DTWireId id(detid);
0123   if (superlayer_to_hits_.find(id.superlayerId().rawId()) == superlayer_to_hits_.end())
0124     return no_hits_;
0125   return superlayer_to_hits_.at(id.superlayerId().rawId());
0126 }
0127
0128 const edm::PSimHitContainer& DTSimHitMatcher::hitsInChamber(unsigned int detid) const {
0129   if (!MuonHitHelper::isDT(detid))
0130     return no_hits_;
0131
0132   const DTWireId id(detid);
0133   if (chamber_to_hits_.find(id.chamberId().rawId()) == chamber_to_hits_.end())
0134     return no_hits_;
0135   return chamber_to_hits_.at(id.chamberId().rawId());
0136 }
0137
0138 bool DTSimHitMatcher::hitStation(int st, int nsuperlayers, int nlayers) const {
0139   int nst = 0;
0140   for (const auto& ddt : chamberIds()) {
0141     const DTChamberId id(ddt);
0142     if (id.station() != st)
0143       continue;
0144
0145     // require at least 1 superlayer
0146     const int nsl(nSuperLayersWithHitsInChamber(id.rawId()));
0147     if (nsl < nsuperlayers)
0148       continue;
0149
0150     // require at least 3 layers hit per chamber
0151     const int nl(nLayersWithHitsInChamber(id.rawId()));
0152     if (nl < nlayers)
0153       continue;
0154     ++nst;
0155   }
0156   return nst;
0157 }
0158
0159 int DTSimHitMatcher::nStations(int nsuperlayers, int nlayers) const {
0160   return (hitStation(1, nsuperlayers, nlayers) + hitStation(2, nsuperlayers, nlayers) +
0161           hitStation(3, nsuperlayers, nlayers) + hitStation(4, nsuperlayers, nlayers));
0162 }
0163
0164 int DTSimHitMatcher::nCellsWithHitsInLayer(unsigned int detid) const {
0165   set<int> layers_with_hits;
0166   const auto& hits = hitsInLayer(detid);
0167   for (const auto& h : hits) {
0168     if (MuonHitHelper::isDT(detid)) {
0169       const DTWireId idd(h.detUnitId());
0170       layers_with_hits.insert(idd.wire());
0171     }
0172   }
0173   return layers_with_hits.size();
0174 }
0175
0176 int DTSimHitMatcher::nLayersWithHitsInSuperLayer(unsigned int detid) const {
0177   set<int> layers_with_hits;
0178   const auto& hits = hitsInSuperLayer(detid);
0179   for (const auto& h : hits) {
0180     if (MuonHitHelper::isDT(detid)) {
0181       const DTLayerId idd(h.detUnitId());
0182       layers_with_hits.insert(idd.layer());
0183     }
0184   }
0185   return layers_with_hits.size();
0186 }
0187
0188 int DTSimHitMatcher::nSuperLayersWithHitsInChamber(unsigned int detid) const {
0189   set<int> sl_with_hits;
0190   const auto& hits = MuonSimHitMatcher::hitsInChamber(detid);
0191   for (const auto& h : hits) {
0192     if (MuonHitHelper::isDT(detid)) {
0193       const DTSuperLayerId idd(h.detUnitId());
0194       sl_with_hits.insert(idd.superLayer());
0195     }
0196   }
0197   return sl_with_hits.size();
0198 }
0199
0200 int DTSimHitMatcher::nLayersWithHitsInChamber(unsigned int detid) const {
0201   int nLayers = 0;
0202   const auto& superLayers(dynamic_cast<const DTGeometry*>(geometry_)->chamber(DTChamberId(detid))->superLayers());
0203   for (const auto& sl : superLayers) {
0204     nLayers += nLayersWithHitsInSuperLayer(sl->id().rawId());
0205   }
0206   return nLayers;
0207 }
0208 float DTSimHitMatcher::simHitsMeanWire(const edm::PSimHitContainer& sim_hits) const {
0209   if (sim_hits.empty())
0210     return -1.f;
0211
0212   float sums = 0.f;
0213   size_t n = 0;
0214   for (const auto& h : sim_hits) {
0215     const LocalPoint& lp = h.entryPoint();
0216     float s;
0217     const auto& d = h.detUnitId();
0218     if (MuonHitHelper::isDT(d)) {
0219       // find nearest wire
0220       s = dynamic_cast<const DTGeometry*>(geometry_)->layer(DTLayerId(d))->specificTopology().channel(lp);
0221     } else
0222       continue;
0223     sums += s;
0224     ++n;
0225   }
0226   if (n == 0)
0227     return -1.f;
0228   return sums / n;
0229 }
0230
0231 std::set<unsigned int> DTSimHitMatcher::hitWiresInDTLayerId(unsigned int detid, int margin_n_wires) const {
0232   set<unsigned int> result;
0233   if (MuonHitHelper::isDT(detid)) {
0234     DTLayerId id(detid);
0235     int max_nwires = dynamic_cast<const DTGeometry*>(geometry_)->layer(id)->specificTopology().channels();
0236     for (int wn = 0; wn <= max_nwires; ++wn) {
0237       DTWireId wid(id, wn);
0238       for (const auto& h : MuonSimHitMatcher::hitsInDetId(wid.rawId())) {
0239         if (verbose_)
0240           edm::LogInfo("DTSimHitMatcher") << "central DTWireId " << wid << " simhit " << h;
0241         int smin = wn - margin_n_wires;
0242         smin = (smin > 0) ? smin : 1;
0243         int smax = wn + margin_n_wires;
0244         smax = (smax <= max_nwires) ? smax : max_nwires;
0245         for (int ss = smin; ss <= smax; ++ss) {
0246           DTWireId widd(id, ss);
0247           if (verbose_)
0248             edm::LogInfo("DTSimHitMatcher") << "\tadding DTWireId to collection " << widd;
0249           result.insert(widd.rawId());
0250         }
0251       }
0252     }
0253   }
0254   return result;
0255 }
0256
0257 std::set<unsigned int> DTSimHitMatcher::hitWiresInDTSuperLayerId(unsigned int detid, int margin_n_wires) const {
0258   set<unsigned int> result;
0259   const auto& layers(dynamic_cast<const DTGeometry*>(geometry_)->superLayer(DTSuperLayerId(detid))->layers());
0260   for (const auto& l : layers) {
0261     if (verbose_)
0262       edm::LogInfo("DTSimHitMatcher") << "hitWiresInDTSuperLayerId::l id " << l->id();
0263     const auto& p(hitWiresInDTLayerId(l->id().rawId(), margin_n_wires));
0264     result.insert(p.begin(), p.end());
0265   }
0266   return result;
0267 }
0268
0269 std::set<unsigned int> DTSimHitMatcher::hitWiresInDTChamberId(unsigned int detid, int margin_n_wires) const {
0270   set<unsigned int> result;
0271   const auto& superLayers(dynamic_cast<const DTGeometry*>(geometry_)->chamber(DTChamberId(detid))->superLayers());
0272   for (const auto& sl : superLayers) {
0273     if (verbose_)
0274       edm::LogInfo("DTSimHitMatcher") << "hitWiresInDTChamberId::sl id " << sl->id();
0275     const auto& p(hitWiresInDTSuperLayerId(sl->id().rawId(), margin_n_wires));
0276     result.insert(p.begin(), p.end());
0277   }
0278   return result;
0279 }
0280
0281 void DTSimHitMatcher::dtChamberIdsToString(const std::set<unsigned int>& set) const {
0282   for (const auto& p : set) {
0283     DTChamberId detId(p);
0284     edm::LogInfo("DTSimHitMatcher") << " " << detId << "\n";
0285   }
0286 }
0287
0288 std::set<unsigned int> DTSimHitMatcher::chamberIdsStation(int station) const {
0289   set<unsigned int> result;
0290   switch (station) {
0291     case 1: {
0292       const auto& p1(chamberIds(MuonHitHelper::DT_MB21p));
0293       const auto& p2(chamberIds(MuonHitHelper::DT_MB11p));
0294       const auto& p3(chamberIds(MuonHitHelper::DT_MB01));
0295       const auto& p4(chamberIds(MuonHitHelper::DT_MB11n));
0296       const auto& p5(chamberIds(MuonHitHelper::DT_MB21n));
0297       result.insert(p1.begin(), p1.end());
0298       result.insert(p2.begin(), p2.end());
0299       result.insert(p3.begin(), p3.end());
0300       result.insert(p4.begin(), p4.end());
0301       result.insert(p5.begin(), p5.end());
0302       break;
0303     }
0304     case 2: {
0305       const auto& p1(chamberIds(MuonHitHelper::DT_MB22p));
0306       const auto& p2(chamberIds(MuonHitHelper::DT_MB12p));
0307       const auto& p3(chamberIds(MuonHitHelper::DT_MB02));
0308       const auto& p4(chamberIds(MuonHitHelper::DT_MB12n));
0309       const auto& p5(chamberIds(MuonHitHelper::DT_MB22n));
0310       result.insert(p1.begin(), p1.end());
0311       result.insert(p2.begin(), p2.end());
0312       result.insert(p3.begin(), p3.end());
0313       result.insert(p4.begin(), p4.end());
0314       result.insert(p5.begin(), p5.end());
0315       break;
0316     }
0317     case 3: {
0318       const auto& p1(chamberIds(MuonHitHelper::DT_MB23p));
0319       const auto& p2(chamberIds(MuonHitHelper::DT_MB13p));
0320       const auto& p3(chamberIds(MuonHitHelper::DT_MB03));
0321       const auto& p4(chamberIds(MuonHitHelper::DT_MB13n));
0322       const auto& p5(chamberIds(MuonHitHelper::DT_MB23n));
0323       result.insert(p1.begin(), p1.end());
0324       result.insert(p2.begin(), p2.end());
0325       result.insert(p3.begin(), p3.end());
0326       result.insert(p4.begin(), p4.end());
0327       result.insert(p5.begin(), p5.end());
0328       break;
0329     }
0330     case 4: {
0331       const auto& p1(chamberIds(MuonHitHelper::DT_MB24p));
0332       const auto& p2(chamberIds(MuonHitHelper::DT_MB14p));
0333       const auto& p3(chamberIds(MuonHitHelper::DT_MB04));
0334       const auto& p4(chamberIds(MuonHitHelper::DT_MB14n));
0335       const auto& p5(chamberIds(MuonHitHelper::DT_MB24n));
0336       result.insert(p1.begin(), p1.end());
0337       result.insert(p2.begin(), p2.end());
0338       result.insert(p3.begin(), p3.end());
0339       result.insert(p4.begin(), p4.end());
0340       result.insert(p5.begin(), p5.end());
0341       break;
0342     }
0343   };
0344   return result;
0345 }