Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:25:26 

0001 #include "SimMuon/MCTruth/interface/MuonAssociatorByHits.h"
0002 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
0003 #include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
0004 #include "DataFormats/DetId/interface/DetId.h"
0005 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0006 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0007 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0008 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0011 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0012 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0013 #include "SimMuon/MCTruth/interface/TrackerMuonHitExtractor.h"
0014 #include <memory>
0015 
0016 #include <sstream>
0017 
0018 using namespace reco;
0019 using namespace std;
0020 using namespace muonAssociatorByHitsDiagnostics;
0021 
0022 namespace muonAssociatorByHitsDiagnostics {
0023   using TrackHitsCollection = MuonAssociatorByHitsHelper::TrackHitsCollection;
0024 
0025   class InputDumper {
0026   public:
0027     InputDumper(const edm::ParameterSet &conf)
0028         : simtracksTag(conf.getParameter<edm::InputTag>("simtracksTag")),
0029           simtracksXFTag(conf.getParameter<edm::InputTag>("simtracksXFTag")),
0030           crossingframe(conf.getParameter<bool>("crossingframe")) {}
0031 
0032     InputDumper(const edm::ParameterSet &conf, edm::ConsumesCollector &&iC) : InputDumper(conf) {
0033       if (crossingframe) {
0034         simtracksXFToken_ = iC.consumes<CrossingFrame<SimTrack>>(simtracksXFTag);
0035         simvertexXFToken_ = iC.consumes<CrossingFrame<SimVertex>>(simtracksXFTag);
0036       } else {
0037         simtracksToken_ = iC.consumes<edm::SimTrackContainer>(simtracksTag);
0038         simvertexToken_ = iC.consumes<edm::SimVertexContainer>(simtracksTag);
0039       }
0040     }
0041 
0042     void dump(const TrackHitsCollection &, const TrackingParticleCollection &, const edm::Event &) const;
0043 
0044   private:
0045     edm::InputTag const simtracksTag;
0046     edm::InputTag const simtracksXFTag;
0047     bool const crossingframe;
0048     edm::EDGetTokenT<CrossingFrame<SimTrack>> simtracksXFToken_;
0049     edm::EDGetTokenT<CrossingFrame<SimVertex>> simvertexXFToken_;
0050     edm::EDGetTokenT<edm::SimTrackContainer> simtracksToken_;
0051     edm::EDGetTokenT<edm::SimVertexContainer> simvertexToken_;
0052   };
0053 
0054   void InputDumper::dump(const TrackHitsCollection &tC,
0055                          const TrackingParticleCollection &tPC,
0056                          const edm::Event &event) const {
0057     // reco::Track collection
0058     edm::LogVerbatim("MuonAssociatorByHits") << "\n"
0059                                              << "reco::Track collection --- size = " << tC.size();
0060 
0061     // TrackingParticle collection
0062     edm::LogVerbatim("MuonAssociatorByHits") << "\n"
0063                                              << "TrackingParticle collection --- size = " << tPC.size();
0064     int j = 0;
0065     for (TrackingParticleCollection::const_iterator ITER = tPC.begin(); ITER != tPC.end(); ITER++, j++) {
0066       edm::LogVerbatim("MuonAssociatorByHits")
0067           << "TrackingParticle " << j << ", q = " << ITER->charge() << ", p = " << ITER->p() << ", pT = " << ITER->pt()
0068           << ", eta = " << ITER->eta() << ", phi = " << ITER->phi();
0069 
0070       edm::LogVerbatim("MuonAssociatorByHits")
0071           << "\t pdg code = " << ITER->pdgId() << ", made of " << ITER->numberOfHits() << " PSimHit"
0072           << " (in " << ITER->numberOfTrackerLayers() << " layers)"
0073           << " from " << ITER->g4Tracks().size() << " SimTrack:";
0074       for (TrackingParticle::g4t_iterator g4T = ITER->g4Track_begin(); g4T != ITER->g4Track_end(); g4T++) {
0075         edm::LogVerbatim("MuonAssociatorByHits") << "\t\t Id:" << g4T->trackId() << "/Evt:(" << g4T->eventId().event()
0076                                                  << "," << g4T->eventId().bunchCrossing() << ")";
0077       }
0078     }
0079 
0080     if (crossingframe) {
0081       const auto &cf_simtracks = event.getHandle(simtracksXFToken_);
0082       unique_ptr<MixCollection<SimTrack>> SimTk(new MixCollection<SimTrack>(cf_simtracks.product()));
0083       edm::LogVerbatim("MuonAssociatorByHits")
0084           << "\n"
0085           << "CrossingFrame<SimTrack> collection with InputTag = " << simtracksXFTag << " has size = " << SimTk->size();
0086       int k = 0;
0087       for (MixCollection<SimTrack>::MixItr ITER = SimTk->begin(); ITER != SimTk->end(); ITER++, k++) {
0088         edm::LogVerbatim("MuonAssociatorByHits")
0089             << "SimTrack " << k << " - Id:" << ITER->trackId() << "/Evt:(" << ITER->eventId().event() << ","
0090             << ITER->eventId().bunchCrossing() << ")"
0091             << " pdgId = " << ITER->type() << ", q = " << ITER->charge() << ", p = " << ITER->momentum().P()
0092             << ", pT = " << ITER->momentum().Pt() << ", eta = " << ITER->momentum().Eta()
0093             << ", phi = " << ITER->momentum().Phi() << "\n * " << *ITER << endl;
0094       }
0095       const auto &cf_simvertices = event.getHandle(simvertexXFToken_);
0096       unique_ptr<MixCollection<SimVertex>> SimVtx(new MixCollection<SimVertex>(cf_simvertices.product()));
0097       edm::LogVerbatim("MuonAssociatorByHits")
0098           << "\n"
0099           << "CrossingFrame<SimVertex> collection with InputTag = " << simtracksXFTag
0100           << " has size = " << SimVtx->size();
0101       int kv = 0;
0102       for (MixCollection<SimVertex>::MixItr VITER = SimVtx->begin(); VITER != SimVtx->end(); VITER++, kv++) {
0103         edm::LogVerbatim("MuonAssociatorByHits") << "SimVertex " << kv << " : " << *VITER << endl;
0104       }
0105     } else {
0106       const auto &simTrackCollection = event.getHandle(simtracksToken_);
0107       const edm::SimTrackContainer simTC = *(simTrackCollection.product());
0108       edm::LogVerbatim("MuonAssociatorByHits")
0109           << "\n"
0110           << "SimTrack collection with InputTag = " << simtracksTag << " has size = " << simTC.size() << endl;
0111       int k = 0;
0112       for (edm::SimTrackContainer::const_iterator ITER = simTC.begin(); ITER != simTC.end(); ITER++, k++) {
0113         edm::LogVerbatim("MuonAssociatorByHits")
0114             << "SimTrack " << k << " - Id:" << ITER->trackId() << "/Evt:(" << ITER->eventId().event() << ","
0115             << ITER->eventId().bunchCrossing() << ")"
0116             << " pdgId = " << ITER->type() << ", q = " << ITER->charge() << ", p = " << ITER->momentum().P()
0117             << ", pT = " << ITER->momentum().Pt() << ", eta = " << ITER->momentum().Eta()
0118             << ", phi = " << ITER->momentum().Phi() << "\n * " << *ITER << endl;
0119       }
0120       const auto &simVertexCollection = event.getHandle(simvertexToken_);
0121       const edm::SimVertexContainer simVC = *(simVertexCollection.product());
0122       edm::LogVerbatim("MuonAssociatorByHits") << "\n"
0123                                                << "SimVertex collection with InputTag = "
0124                                                << "g4SimHits"
0125                                                << " has size = " << simVC.size() << endl;
0126       int kv = 0;
0127       for (edm::SimVertexContainer::const_iterator VITER = simVC.begin(); VITER != simVC.end(); VITER++, kv++) {
0128         edm::LogVerbatim("MuonAssociatorByHits") << "SimVertex " << kv << " : " << *VITER << endl;
0129       }
0130     }
0131   }
0132 
0133 }  // namespace muonAssociatorByHitsDiagnostics
0134 
0135 MuonAssociatorByHits::MuonAssociatorByHits(const edm::ParameterSet &conf, edm::ConsumesCollector &&iC)
0136     : helper_(conf),
0137       trackerHitAssociatorConfig_(conf, std::move(iC)),
0138       gemHitAssociatorConfig_(conf, iC),
0139       rpcHitAssociatorConfig_(conf, iC),
0140       cscHitAssociatorConfig_(conf, iC),
0141       dtHitAssociatorConfig_(conf, iC),
0142       ttopoToken_(iC.esConsumes()) {
0143   // hack for consumes
0144   if (conf.getUntrackedParameter<bool>("dumpInputCollections")) {
0145     diagnostics_ = std::make_unique<InputDumper>(conf, std::move(iC));
0146   }
0147 }
0148 
0149 MuonAssociatorByHits::~MuonAssociatorByHits() {}
0150 
0151 RecoToSimCollection MuonAssociatorByHits::associateRecoToSim(
0152     const edm::RefToBaseVector<reco::Track> &tC,
0153     const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
0154     const edm::Event *e,
0155     const edm::EventSetup *setup) const {
0156   RecoToSimCollection outputCollection(&e->productGetter());
0157 
0158   MuonAssociatorByHitsHelper::TrackHitsCollection tH;
0159   for (auto it = tC.begin(), ed = tC.end(); it != ed; ++it) {
0160     tH.push_back(std::make_pair((*it)->recHitsBegin(), (*it)->recHitsEnd()));
0161   }
0162 
0163   // Retrieve tracker topology from geometry
0164   const TrackerTopology *tTopo = &setup->getData(ttopoToken_);
0165 
0166   // Tracker hit association
0167   TrackerHitAssociator trackertruth(*e, trackerHitAssociatorConfig_);
0168   // CSC hit association
0169   CSCHitAssociator csctruth(*e, *setup, cscHitAssociatorConfig_);
0170   // DT hit association
0171   bool printRtS(true);
0172   DTHitAssociator dttruth(*e, *setup, dtHitAssociatorConfig_, printRtS);
0173   // RPC hit association
0174   RPCHitAssociator rpctruth(*e, rpcHitAssociatorConfig_);
0175   // GEM hit association
0176   GEMHitAssociator gemtruth(*e, gemHitAssociatorConfig_);
0177 
0178   MuonAssociatorByHitsHelper::Resources resources = {
0179       tTopo, &trackertruth, &csctruth, &dttruth, &rpctruth, &gemtruth, {}};
0180 
0181   if (diagnostics_) {
0182     resources.diagnostics_ = [this, e](const TrackHitsCollection &hC, const TrackingParticleCollection &pC) {
0183       diagnostics_->dump(hC, pC, *e);
0184     };
0185   }
0186 
0187   auto bareAssoc = helper_.associateRecoToSimIndices(tH, TPCollectionH, resources);
0188   for (auto it = bareAssoc.begin(), ed = bareAssoc.end(); it != ed; ++it) {
0189     for (auto itma = it->second.begin(), edma = it->second.end(); itma != edma; ++itma) {
0190       outputCollection.insert(tC[it->first], std::make_pair(TPCollectionH[itma->idx], itma->quality));
0191     }
0192   }
0193 
0194   outputCollection.post_insert();  // perhaps not even necessary
0195   return outputCollection;
0196 }
0197 
0198 SimToRecoCollection MuonAssociatorByHits::associateSimToReco(
0199     const edm::RefToBaseVector<reco::Track> &tC,
0200     const edm::RefVector<TrackingParticleCollection> &TPCollectionH,
0201     const edm::Event *e,
0202     const edm::EventSetup *setup) const {
0203   SimToRecoCollection outputCollection(&e->productGetter());
0204   MuonAssociatorByHitsHelper::TrackHitsCollection tH;
0205   for (auto it = tC.begin(), ed = tC.end(); it != ed; ++it) {
0206     tH.push_back(std::make_pair((*it)->recHitsBegin(), (*it)->recHitsEnd()));
0207   }
0208 
0209   // Retrieve tracker topology from geometry
0210   const TrackerTopology *tTopo = &setup->getData(ttopoToken_);
0211 
0212   // Tracker hit association
0213   TrackerHitAssociator trackertruth(*e, trackerHitAssociatorConfig_);
0214   // CSC hit association
0215   CSCHitAssociator csctruth(*e, *setup, cscHitAssociatorConfig_);
0216   // DT hit association
0217   bool printRtS = false;
0218   DTHitAssociator dttruth(*e, *setup, dtHitAssociatorConfig_, printRtS);
0219   // RPC hit association
0220   RPCHitAssociator rpctruth(*e, rpcHitAssociatorConfig_);
0221   // GEM hit association
0222   GEMHitAssociator gemtruth(*e, gemHitAssociatorConfig_);
0223 
0224   MuonAssociatorByHitsHelper::Resources resources = {
0225       tTopo, &trackertruth, &csctruth, &dttruth, &rpctruth, &gemtruth, {}};
0226 
0227   auto bareAssoc = helper_.associateSimToRecoIndices(tH, TPCollectionH, resources);
0228   for (auto it = bareAssoc.begin(), ed = bareAssoc.end(); it != ed; ++it) {
0229     for (auto itma = it->second.begin(), edma = it->second.end(); itma != edma; ++itma) {
0230       outputCollection.insert(TPCollectionH[it->first], std::make_pair(tC[itma->idx], itma->quality));
0231     }
0232   }
0233 
0234   outputCollection.post_insert();  // perhaps not even necessary
0235   return outputCollection;
0236 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!