MuonHits/src/CSCSimHitMatcher.cc

0001 #include "Validation/MuonHits/interface/CSCSimHitMatcher.h"
0002 #include "TGraphErrors.h"
0003 #include "TF1.h"
0004 #include "TMinuitMinimizer.h"
0005
0006 using namespace std;
0007
0008 CSCSimHitMatcher::CSCSimHitMatcher(const edm::ParameterSet& ps, edm::ConsumesCollector&& iC)
0009     : MuonSimHitMatcher(ps, std::move(iC)) {
0010   simHitPSet_ = ps.getParameterSet("cscSimHit");
0011   verbose_ = simHitPSet_.getParameter<int>("verbose");
0012   simMuOnly_ = simHitPSet_.getParameter<bool>("simMuOnly");
0013   discardEleHits_ = simHitPSet_.getParameter<bool>("discardEleHits");
0014
0015   simHitInput_ = iC.consumes<edm::PSimHitContainer>(simHitPSet_.getParameter<edm::InputTag>("inputTag"));
0016   geomToken_ = iC.esConsumes<CSCGeometry, MuonGeometryRecord>();
0017
0018   //In order to make fitting ROOT histograms thread safe
0019   // one must call this undocumented function
0020   TMinuitMinimizer::UseStaticMinuit(false);
0021 }
0022
0023 /// initialize the event
0024 void CSCSimHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0025   geometry_ = &iSetup.getData(geomToken_);
0026   MuonSimHitMatcher::init(iEvent, iSetup);
0027 }
0028
0029 /// do the matching
0030 void CSCSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) {
0031   clear();
0032
0033   // instantiates the track ids and simhits
0034   MuonSimHitMatcher::match(track, vertex);
0035
0036   // hard cut on non-CSC muons
0037   if (std::abs(track.momentum().eta()) < 0.9)
0038     return;
0039   if (std::abs(track.momentum().eta()) > 2.45)
0040     return;
0041
0042   matchSimHitsToSimTrack();
0043
0044   if (verbose_) {
0045     edm::LogInfo("CSCSimHitMatcher") << "nTrackIds " << track_ids_.size() << " nSelectedCSCSimHits " << hits_.size();
0046     edm::LogInfo("CSCSimHitMatcher") << "detids CSC " << detIds(0).size();
0047
0048     for (const auto& id : detIds(0)) {
0049       const auto& simhits = hitsInDetId(id);
0050       const auto& simhits_gp = simHitsMeanPosition(simhits);
0051       const auto& strips = hitStripsInDetId(id);
0052       CSCDetId cscid(id);
0053       if (cscid.station() == 1 and (cscid.ring() == 1 or cscid.ring() == 4)) {
0054         edm::LogInfo("CSCSimHitMatcher") << "cscdetid " << CSCDetId(id) << ": " << simhits.size() << " "
0055                                          << simhits_gp.phi() << " " << detid_to_hits_[id].size();
0056         edm::LogInfo("CSCSimHitMatcher") << "nStrip " << strips.size();
0057         edm::LogInfo("CSCSimHitMatcher") << "strips : ";
0058         for (const auto& p : strips) {
0059           edm::LogInfo("CSCSimHitMatcher") << p;
0060         }
0061       }
0062     }
0063   }
0064 }
0065
0066 void CSCSimHitMatcher::matchSimHitsToSimTrack() {
0067   for (const auto& track_id : track_ids_) {
0068     for (const auto& h : simHits_) {
0069       if (h.trackId() != track_id)
0070         continue;
0071       int pdgid = h.particleType();
0072       if (simMuOnly_ && std::abs(pdgid) != 13)
0073         continue;
0074       // discard electron hits in the CSC chambers
0075       if (discardEleHits_ && std::abs(pdgid) == 11)
0076         continue;
0077
0078       const CSCDetId& layer_id(h.detUnitId());
0079       hits_.push_back(h);
0080       detid_to_hits_[h.detUnitId()].push_back(h);
0081       chamber_to_hits_[layer_id.chamberId().rawId()].push_back(h);
0082     }
0083   }
0084 }
0085
0086 std::set<unsigned int> CSCSimHitMatcher::detIds(int csc_type) const {
0087   std::set<unsigned int> result;
0088   for (const auto& p : detid_to_hits_) {
0089     const auto& id = p.first;
0090     if (csc_type > 0) {
0091       CSCDetId detId(id);
0092       if (MuonHitHelper::toCSCType(detId.station(), detId.ring()) != csc_type)
0093         continue;
0094     }
0095     result.insert(id);
0096   }
0097   return result;
0098 }
0099
0100 std::set<unsigned int> CSCSimHitMatcher::chamberIds(int csc_type) const {
0101   std::set<unsigned int> result;
0102   for (const auto& p : chamber_to_hits_) {
0103     const auto& id = p.first;
0104     if (csc_type > 0) {
0105       CSCDetId detId(id);
0106       if (MuonHitHelper::toCSCType(detId.station(), detId.ring()) != csc_type)
0107         continue;
0108     }
0109     result.insert(id);
0110   }
0111   return result;
0112 }
0113
0114 int CSCSimHitMatcher::nLayersWithHitsInChamber(unsigned int detid) const {
0115   set<int> layers_with_hits;
0116   for (const auto& h : MuonSimHitMatcher::hitsInChamber(detid)) {
0117     const CSCDetId& idd(h.detUnitId());
0118     layers_with_hits.insert(idd.layer());
0119   }
0120   return layers_with_hits.size();
0121 }
0122
0123 bool CSCSimHitMatcher::hitStation(int st, int nlayers) const {
0124   int nst = 0;
0125   for (const auto& ddt : chamberIds()) {
0126     const CSCDetId id(ddt);
0127     int ri(id.ring());
0128     if (id.station() != st)
0129       continue;
0130
0131     // ME1/a case - check the ME1/b chamber
0132     if (st == 1 and ri == 4) {
0133       CSCDetId idME1b(id.endcap(), id.station(), 1, id.chamber(), 0);
0134       const int nl1a(nLayersWithHitsInChamber(id.rawId()));
0135       const int nl1b(nLayersWithHitsInChamber(idME1b.rawId()));
0136       if (nl1a + nl1b < nlayers)
0137         continue;
0138       ++nst;
0139     }
0140     // ME1/b case - check the ME1/a chamber
0141     else if (st == 1 and ri == 1) {
0142       CSCDetId idME1a(id.endcap(), id.station(), 4, id.chamber(), 0);
0143       const int nl1a(nLayersWithHitsInChamber(idME1a.rawId()));
0144       const int nl1b(nLayersWithHitsInChamber(id.rawId()));
0145       if (nl1a + nl1b < nlayers)
0146         continue;
0147       ++nst;
0148     }
0149     // default case
0150     else {
0151       const int nl(nLayersWithHitsInChamber(id.rawId()));
0152       if (nl < nlayers)
0153         continue;
0154       ++nst;
0155     }
0156   }
0157   return nst;
0158 }
0159
0160 int CSCSimHitMatcher::nStations(int nlayers) const {
0161   return (hitStation(1, nlayers) + hitStation(2, nlayers) + hitStation(3, nlayers) + hitStation(4, nlayers));
0162 }
0163
0164 float CSCSimHitMatcher::LocalBendingInChamber(unsigned int detid) const {
0165   const CSCDetId cscid(detid);
0166   if (nLayersWithHitsInChamber(detid) < 6)
0167     return -100;
0168   float phi_layer1 = -10;
0169   float phi_layer6 = 10;
0170
0171   if (cscid.station() == 1 and (cscid.ring() == 1 or cscid.ring() == 4)) {
0172     // phi in layer 1
0173     const CSCDetId cscid1a(cscid.endcap(), cscid.station(), 4, cscid.chamber(), 1);
0174     const CSCDetId cscid1b(cscid.endcap(), cscid.station(), 1, cscid.chamber(), 1);
0175     const edm::PSimHitContainer& hits1a = MuonSimHitMatcher::hitsInDetId(cscid1a.rawId());
0176     const edm::PSimHitContainer& hits1b = MuonSimHitMatcher::hitsInDetId(cscid1b.rawId());
0177     const GlobalPoint& gp1a = simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid1a.rawId()));
0178     const GlobalPoint& gp1b = simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid1b.rawId()));
0179     if (!hits1a.empty() and !hits1b.empty())
0180       phi_layer1 = (gp1a.phi() + gp1b.phi()) / 2.0;
0181     else if (!hits1a.empty())
0182       phi_layer1 = gp1a.phi();
0183     else if (!hits1b.empty())
0184       phi_layer1 = gp1b.phi();
0185
0186     // phi in layer 6
0187     const CSCDetId cscid6a(cscid.endcap(), cscid.station(), 4, cscid.chamber(), 6);
0188     const CSCDetId cscid6b(cscid.endcap(), cscid.station(), 1, cscid.chamber(), 6);
0189     const edm::PSimHitContainer& hits6a = MuonSimHitMatcher::hitsInDetId(cscid6a.rawId());
0190     const edm::PSimHitContainer& hits6b = MuonSimHitMatcher::hitsInDetId(cscid6b.rawId());
0191     const GlobalPoint& gp6a = simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid6a.rawId()));
0192     const GlobalPoint& gp6b = simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid6b.rawId()));
0193     if (!hits6a.empty() and !hits6b.empty())
0194       phi_layer6 = (gp6a.phi() + gp6b.phi()) / 2.0;
0195     else if (!hits6a.empty())
0196       phi_layer6 = gp6a.phi();
0197     else if (!hits6b.empty())
0198       phi_layer6 = gp6b.phi();
0199
0200   } else {
0201     // phi in layer 1
0202     const CSCDetId cscid1(cscid.endcap(), cscid.station(), cscid.ring(), cscid.chamber(), 1);
0203     const GlobalPoint& gp1 = simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid1.rawId()));
0204     phi_layer1 = gp1.phi();
0205
0206     // phi in layer 6
0207     const CSCDetId cscid6(cscid.endcap(), cscid.station(), cscid.ring(), cscid.chamber(), 6);
0208     const GlobalPoint& gp6 = simHitsMeanPosition(MuonSimHitMatcher::hitsInDetId(cscid6.rawId()));
0209     phi_layer6 = gp6.phi();
0210   }
0211   return deltaPhi(phi_layer6, phi_layer1);
0212 }
0213
0214 // difference in strip per layer
0215 void CSCSimHitMatcher::fitHitsInChamber(unsigned int detid, float& intercept, float& slope) const {
0216   const CSCDetId cscid(detid);
0217
0218   const auto& sim_hits = hitsInChamber(detid);
0219
0220   if (sim_hits.empty())
0221     return;
0222
0223   vector<float> x;
0224   vector<float> y;
0225   vector<float> xe;
0226   vector<float> ye;
0227
0228   const float HALF_STRIP_ERROR = 0.288675;
0229
0230   for (const auto& h : sim_hits) {
0231     const LocalPoint& lp = h.entryPoint();
0232     const auto& d = h.detUnitId();
0233     float s = dynamic_cast<const CSCGeometry*>(geometry_)->layer(d)->geometry()->strip(lp);
0234     // shift to key half strip layer (layer 3)
0235     x.push_back(CSCDetId(d).layer() - 3);
0236     y.push_back(s);
0237     xe.push_back(float(0));
0238     ye.push_back(2 * HALF_STRIP_ERROR);
0239   }
0240   if (x.size() < 2)
0241     return;
0242
0243   std::unique_ptr<TGraphErrors> gr(new TGraphErrors(x.size(), &x[0], &y[0], &xe[0], &ye[0]));
0244   TF1 fit("fit", "pol1", -3, 4);
0245   gr->Fit(&fit, "EMQN");
0246
0247   intercept = fit.GetParameter(0);
0248   slope = fit.GetParameter(1);
0249 }
0250
0251 float CSCSimHitMatcher::simHitsMeanStrip(const edm::PSimHitContainer& sim_hits) const {
0252   if (sim_hits.empty())
0253     return -1.f;
0254
0255   float sums = 0.f;
0256   size_t n = 0;
0257   for (const auto& h : sim_hits) {
0258     const LocalPoint& lp = h.entryPoint();
0259     float s;
0260     const auto& d = h.detUnitId();
0261     s = dynamic_cast<const CSCGeometry*>(geometry_)->layer(d)->geometry()->strip(lp);
0262     // convert to half-strip:
0263     s *= 2.;
0264     sums += s;
0265     ++n;
0266   }
0267   if (n == 0)
0268     return -1.f;
0269   return sums / n;
0270 }
0271
0272 float CSCSimHitMatcher::simHitsMeanWG(const edm::PSimHitContainer& sim_hits) const {
0273   if (sim_hits.empty())
0274     return -1.f;
0275
0276   float sums = 0.f;
0277   size_t n = 0;
0278   for (const auto& h : sim_hits) {
0279     const LocalPoint& lp = h.entryPoint();
0280     float s;
0281     const auto& d = h.detUnitId();
0282     // find nearest wire
0283     const auto& layerG(dynamic_cast<const CSCGeometry*>(geometry_)->layer(d)->geometry());
0284     int nearestWire(layerG->nearestWire(lp));
0285     // then find the corresponding wire group
0286     s = layerG->wireGroup(nearestWire);
0287     sums += s;
0288     ++n;
0289   }
0290   if (n == 0)
0291     return -1.f;
0292   return sums / n;
0293 }
0294
0295 std::set<int> CSCSimHitMatcher::hitStripsInDetId(unsigned int detid, int margin_n_strips) const {
0296   set<int> result;
0297   const auto& simhits = MuonSimHitMatcher::hitsInDetId(detid);
0298   CSCDetId id(detid);
0299   int max_nstrips = dynamic_cast<const CSCGeometry*>(geometry_)->layer(id)->geometry()->numberOfStrips();
0300   for (const auto& h : simhits) {
0301     const LocalPoint& lp = h.entryPoint();
0302     int central_strip = dynamic_cast<const CSCGeometry*>(geometry_)->layer(id)->geometry()->nearestStrip(lp);
0303     int smin = central_strip - margin_n_strips;
0304     smin = (smin > 0) ? smin : 1;
0305     int smax = central_strip + margin_n_strips;
0306     smax = (smax <= max_nstrips) ? smax : max_nstrips;
0307     for (int ss = smin; ss <= smax; ++ss)
0308       result.insert(ss);
0309   }
0310   return result;
0311 }
0312
0313 std::set<int> CSCSimHitMatcher::hitWiregroupsInDetId(unsigned int detid, int margin_n_wg) const {
0314   set<int> result;
0315   const auto& simhits = MuonSimHitMatcher::hitsInDetId(detid);
0316   CSCDetId id(detid);
0317   int max_n_wg = dynamic_cast<const CSCGeometry*>(geometry_)->layer(id)->geometry()->numberOfWireGroups();
0318   for (const auto& h : simhits) {
0319     const LocalPoint& lp = h.entryPoint();
0320     const auto& layer_geo = dynamic_cast<const CSCGeometry*>(geometry_)->layer(id)->geometry();
0321     int central_wg = layer_geo->wireGroup(layer_geo->nearestWire(lp));
0322     int wg_min = central_wg - margin_n_wg;
0323     wg_min = (wg_min > 0) ? wg_min : 1;
0324     int wg_max = central_wg + margin_n_wg;
0325     wg_max = (wg_max <= max_n_wg) ? wg_max : max_n_wg;
0326     for (int wg = wg_min; wg <= wg_max; ++wg)
0327       result.insert(wg);
0328   }
0329   return result;
0330 }
0331
0332 int CSCSimHitMatcher::nCoincidenceChambers(int min_n_layers) const {
0333   int result = 0;
0334   const auto& chamber_ids = chamberIds(0);
0335   for (const auto& id : chamber_ids) {
0336     if (nLayersWithHitsInChamber(id) >= min_n_layers)
0337       result += 1;
0338   }
0339   return result;
0340 }
0341
0342 void CSCSimHitMatcher::chamberIdsToString(const std::set<unsigned int>& set) const {
0343   for (const auto& p : set) {
0344     CSCDetId detId(p);
0345     edm::LogInfo("CSCSimHitMatcher") << " " << detId << "\n";
0346   }
0347 }
0348
0349 std::set<unsigned int> CSCSimHitMatcher::chamberIdsStation(int station) const {
0350   set<unsigned int> result;
0351   switch (station) {
0352     case 1: {
0353       const auto& p1(chamberIds(MuonHitHelper::CSC_ME1a));
0354       const auto& p2(chamberIds(MuonHitHelper::CSC_ME1b));
0355       const auto& p3(chamberIds(MuonHitHelper::CSC_ME12));
0356       const auto& p4(chamberIds(MuonHitHelper::CSC_ME13));
0357       result.insert(p1.begin(), p1.end());
0358       result.insert(p2.begin(), p2.end());
0359       result.insert(p3.begin(), p3.end());
0360       result.insert(p4.begin(), p4.end());
0361       break;
0362     }
0363     case 2: {
0364       const auto& p1(chamberIds(MuonHitHelper::CSC_ME21));
0365       const auto& p2(chamberIds(MuonHitHelper::CSC_ME22));
0366       result.insert(p1.begin(), p1.end());
0367       result.insert(p2.begin(), p2.end());
0368       break;
0369     }
0370     case 3: {
0371       const auto& p1(chamberIds(MuonHitHelper::CSC_ME31));
0372       const auto& p2(chamberIds(MuonHitHelper::CSC_ME32));
0373       result.insert(p1.begin(), p1.end());
0374       result.insert(p2.begin(), p2.end());
0375       break;
0376     }
0377     case 4: {
0378       const auto& p1(chamberIds(MuonHitHelper::CSC_ME41));
0379       const auto& p2(chamberIds(MuonHitHelper::CSC_ME42));
0380       result.insert(p1.begin(), p1.end());
0381       result.insert(p2.begin(), p2.end());
0382       break;
0383     }
0384   };
0385   return result;
0386 }
0387
0388 void CSCSimHitMatcher::clear() { MuonSimHitMatcher::clear(); }