MuonHits/src/MuonSimHitMatcher.cc

0001 #include "Validation/MuonHits/interface/MuonSimHitMatcher.h"
0002
0003 #include <algorithm>
0004
0005 using namespace std;
0006
0007 MuonSimHitMatcher::MuonSimHitMatcher(const edm::ParameterSet& ps, edm::ConsumesCollector&& iC) {
0008   const auto& simVertex = ps.getParameterSet("simVertex");
0009   const auto& simTrack = ps.getParameterSet("simTrack");
0010   verboseSimTrack_ = simTrack.getParameter<int>("verbose");
0011
0012   simVertexInput_ = iC.consumes<edm::SimVertexContainer>(simVertex.getParameter<edm::InputTag>("inputTag"));
0013   simTrackInput_ = iC.consumes<edm::SimTrackContainer>(simTrack.getParameter<edm::InputTag>("inputTag"));
0014 }
0015
0016 /// initialize the event
0017 void MuonSimHitMatcher::init(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0018   hasGeometry_ = true;
0019
0020   iEvent.getByToken(simTrackInput_, simTracksH_);
0021   iEvent.getByToken(simVertexInput_, simVerticesH_);
0022   iEvent.getByToken(simHitInput_, simHitsH_);
0023 }
0024
0025 /// do the matching
0026 void MuonSimHitMatcher::match(const SimTrack& track, const SimVertex& vertex) {
0027   simTracks_ = *simTracksH_.product();
0028   simVertices_ = *simVerticesH_.product();
0029   simHits_ = *simHitsH_.product();
0030
0031   // fill trkId2Index association:
0032   int no = 0;
0033   trkid_to_index_.clear();
0034   for (const auto& t : simTracks_) {
0035     trkid_to_index_[t.trackId()] = no;
0036     no++;
0037   }
0038
0039   track_ids_ = getIdsOfSimTrackShower(track.trackId(), simTracks_, simVertices_);
0040   if (verboseSimTrack_) {
0041     edm::LogInfo("MuonSimHitMatcher") << "Printing track_ids" << std::endl;
0042     for (const auto& id : track_ids_)
0043       edm::LogInfo("MuonSimHitMatcher") << "id: " << id << std::endl;
0044   }
0045 }
0046
0047 std::vector<unsigned int> MuonSimHitMatcher::getIdsOfSimTrackShower(unsigned int initial_trk_id,
0048                                                                     const edm::SimTrackContainer& simTracks,
0049                                                                     const edm::SimVertexContainer& simVertices) {
0050   vector<unsigned int> result;
0051   result.push_back(initial_trk_id);
0052
0053   if (!simMuOnly_)
0054     return result;
0055
0056   for (const auto& t : simTracks_) {
0057     SimTrack last_trk = t;
0058     // if (std::abs(t.type()) != 13) continue;
0059     bool is_child = false;
0060     while (true) {
0061       if (last_trk.noVertex())
0062         break;
0063       if (simVertices_[last_trk.vertIndex()].noParent())
0064         break;
0065
0066       unsigned parentId = simVertices_[last_trk.vertIndex()].parentIndex();
0067       if (parentId == initial_trk_id) {
0068         is_child = true;
0069         break;
0070       }
0071
0072       const auto& association = trkid_to_index_.find(parentId);
0073       if (association == trkid_to_index_.end())
0074         break;
0075
0076       last_trk = simTracks_[association->second];
0077     }
0078     if (is_child) {
0079       result.push_back(t.trackId());
0080     }
0081   }
0082   return result;
0083 }
0084
0085 const edm::PSimHitContainer& MuonSimHitMatcher::simHits(int sub) const { return hits_; }
0086
0087 const edm::PSimHitContainer& MuonSimHitMatcher::hitsInDetId(unsigned int detid) const {
0088   if (detid_to_hits_.find(detid) == detid_to_hits_.end())
0089     return no_hits_;
0090   return detid_to_hits_.at(detid);
0091 }
0092
0093 const edm::PSimHitContainer& MuonSimHitMatcher::hitsInChamber(unsigned int detid) const {
0094   if (chamber_to_hits_.find(detid) == chamber_to_hits_.end())
0095     return no_hits_;
0096   return chamber_to_hits_.at(detid);
0097 }
0098
0099 GlobalPoint MuonSimHitMatcher::simHitsMeanPosition(const edm::PSimHitContainer& sim_hits) const {
0100   if (sim_hits.empty())
0101     return GlobalPoint();  // point "zero"
0102
0103   float sumx, sumy, sumz;
0104   sumx = sumy = sumz = 0.f;
0105   size_t n = 0;
0106   for (const auto& h : sim_hits) {
0107     const LocalPoint& lp = h.entryPoint();
0108     const GlobalPoint& gp = geometry_->idToDet(h.detUnitId())->surface().toGlobal(lp);
0109     sumx += gp.x();
0110     sumy += gp.y();
0111     sumz += gp.z();
0112     ++n;
0113   }
0114   if (n == 0)
0115     return GlobalPoint();
0116   return GlobalPoint(sumx / n, sumy / n, sumz / n);
0117 }
0118
0119 GlobalVector MuonSimHitMatcher::simHitsMeanMomentum(const edm::PSimHitContainer& sim_hits) const {
0120   if (sim_hits.empty())
0121     return GlobalVector();  // point "zero"
0122
0123   float sumx, sumy, sumz;
0124   sumx = sumy = sumz = 0.f;
0125   size_t n = 0;
0126   for (const auto& h : sim_hits) {
0127     const LocalVector& lv = h.momentumAtEntry();
0128     const GlobalVector& gv = geometry_->idToDet(h.detUnitId())->surface().toGlobal(lv);
0129     sumx += gv.x();
0130     sumy += gv.y();
0131     sumz += gv.z();
0132     ++n;
0133   }
0134   if (n == 0)
0135     return GlobalVector();
0136   return GlobalVector(sumx / n, sumy / n, sumz / n);
0137 }
0138
0139 void MuonSimHitMatcher::clear() {
0140   track_ids_.clear();
0141   trkid_to_index_.clear();
0142   hits_.clear();
0143   detid_to_hits_.clear();
0144   chamber_to_hits_.clear();
0145 }