MuonToSimAssociatorByHits.cc

Go to the documentation of this file.

#include "DataFormats/CSCRecHit/interface/CSCSegment.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "SimMuon/MCTruth/interface/MuonToSimAssociatorByHits.h"
#include "SimMuon/MCTruth/interface/TrackerMuonHitExtractor.h"
#include <sstream>

using namespace reco;
using namespace std;

MuonToSimAssociatorByHits::MuonToSimAssociatorByHits(const edm::ParameterSet &conf, edm::ConsumesCollector &&iC)
    : helper_(conf), conf_(conf), trackerHitAssociatorConfig_(conf, std::move(iC)) {
  TrackerMuonHitExtractor hitExtractor(conf_, std::move(iC));

  // hack for consumes
  RPCHitAssociator rpctruth(conf, std::move(iC));
  GEMHitAssociator gemtruth(conf, std::move(iC));
  DTHitAssociator dttruth(conf, std::move(iC));
  CSCHitAssociator muonTruth(conf, std::move(iC));
}

MuonToSimAssociatorByHits::~MuonToSimAssociatorByHits() {}

void MuonToSimAssociatorByHits::associateMuons(MuonToSimCollection &recToSim,
                                               SimToMuonCollection &simToRec,
                                               const edm::Handle<edm::View<reco::Muon>> &tCH,
                                               MuonTrackType type,
                                               const edm::Handle<TrackingParticleCollection> &tPCH,
                                               const edm::Event *event,
                                               const edm::EventSetup *setup) const {
  edm::RefVector<TrackingParticleCollection> tpc(tPCH.id());
  for (unsigned int j = 0; j < tPCH->size(); j++)
    tpc.push_back(edm::Ref<TrackingParticleCollection>(tPCH, j));

  edm::RefToBaseVector<reco::Muon> muonBaseRefVector;
  for (size_t i = 0; i < tCH->size(); ++i)
    muonBaseRefVector.push_back(tCH->refAt(i));

  associateMuons(recToSim, simToRec, muonBaseRefVector, type, tpc, event, setup);
}

void MuonToSimAssociatorByHits::associateMuons(MuonToSimCollection &recToSim,
                                               SimToMuonCollection &simToRec,
                                               const edm::RefToBaseVector<reco::Muon> &muons,
                                               MuonTrackType trackType,
                                               const edm::RefVector<TrackingParticleCollection> &tPC,
                                               const edm::Event *event,
                                               const edm::EventSetup *setup) const {
  MuonAssociatorByHitsHelper::TrackHitsCollection muonHitRefs;
  edm::OwnVector<TrackingRecHit> allTMRecHits;  // this I will fill in only for tracker muon hits from
                                                // segments
  switch (trackType) {
    case InnerTk:
      for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
        edm::RefToBase<reco::Muon> mur = *it;
        if (mur->track().isNonnull()) {
          muonHitRefs.push_back(std::make_pair(mur->track()->recHitsBegin(), mur->track()->recHitsEnd()));
        } else {
          muonHitRefs.push_back(std::make_pair(allTMRecHits.data().end(), allTMRecHits.data().end()));
        }
      }
      break;
    case OuterTk:
      for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
        edm::RefToBase<reco::Muon> mur = *it;
        if (mur->outerTrack().isNonnull()) {
          muonHitRefs.push_back(std::make_pair(mur->outerTrack()->recHitsBegin(), mur->outerTrack()->recHitsEnd()));
        } else {
          muonHitRefs.push_back(std::make_pair(allTMRecHits.data().end(), allTMRecHits.data().end()));
        }
      }
      break;
    case GlobalTk:
      for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
        edm::RefToBase<reco::Muon> mur = *it;
        if (mur->globalTrack().isNonnull()) {
          muonHitRefs.push_back(std::make_pair(mur->globalTrack()->recHitsBegin(), mur->globalTrack()->recHitsEnd()));
        } else {
          muonHitRefs.push_back(std::make_pair(allTMRecHits.data().end(), allTMRecHits.data().end()));
        }
      }
      break;
    case Segments: {
      TrackerMuonHitExtractor hitExtractor(conf_);
      hitExtractor.init(*event, *setup);
      // puts hits in the vector, and record indices
      std::vector<std::pair<size_t, size_t>> muonHitIndices;
      for (edm::RefToBaseVector<reco::Muon>::const_iterator it = muons.begin(), ed = muons.end(); it != ed; ++it) {
        edm::RefToBase<reco::Muon> mur = *it;
        std::pair<size_t, size_t> indices(allTMRecHits.size(), allTMRecHits.size());
        if (mur->isTrackerMuon()) {
          std::vector<const TrackingRecHit *> hits = hitExtractor.getMuonHits(*mur);
          for (std::vector<const TrackingRecHit *>::const_iterator ith = hits.begin(), edh = hits.end(); ith != edh;
               ++ith) {
            allTMRecHits.push_back(**ith);
          }
          indices.second += hits.size();
        }
        muonHitIndices.push_back(indices);
      }
      // convert indices into pairs of iterators to references
      typedef std::pair<size_t, size_t> index_pair;
      trackingRecHit_iterator hitRefBegin = allTMRecHits.data().begin();
      for (std::vector<std::pair<size_t, size_t>>::const_iterator idxs = muonHitIndices.begin(),
                                                                  idxend = muonHitIndices.end();
           idxs != idxend;
           ++idxs) {
        muonHitRefs.push_back(std::make_pair(hitRefBegin + idxs->first, hitRefBegin + idxs->second));
      }

    } break;
  }

  // Retrieve tracker topology from geometry
  edm::ESHandle<TrackerTopology> tTopoHand;
  setup->get<TrackerTopologyRcd>().get(tTopoHand);
  const TrackerTopology *tTopo = tTopoHand.product();

  // Tracker hit association
  TrackerHitAssociator trackertruth(*event, trackerHitAssociatorConfig_);
  // CSC hit association
  CSCHitAssociator csctruth(*event, *setup, conf_);
  // DT hit association
  bool printRtS = true;
  DTHitAssociator dttruth(*event, *setup, conf_, printRtS);
  // RPC hit association
  RPCHitAssociator rpctruth(*event, *setup, conf_);
  // GEM hit association
  GEMHitAssociator gemtruth(*event, *setup, conf_);

  MuonAssociatorByHitsHelper::Resources resources = {tTopo, &trackertruth, &csctruth, &dttruth, &rpctruth, &gemtruth};

  auto recSimColl = helper_.associateRecoToSimIndices(muonHitRefs, tPC, resources);
  for (auto it = recSimColl.begin(), ed = recSimColl.end(); it != ed; ++it) {
    edm::RefToBase<reco::Muon> rec = muons[it->first];
    std::vector<std::pair<TrackingParticleRef, double>> &tpAss = recToSim[rec];
    for (auto const &a : it->second) {
      tpAss.push_back(std::make_pair(tPC[a.idx], a.quality));
    }
  }
  auto simRecColl = helper_.associateSimToRecoIndices(muonHitRefs, tPC, resources);
  for (auto it = simRecColl.begin(), ed = simRecColl.end(); it != ed; ++it) {
    TrackingParticleRef sim = tPC[it->first];
    std::vector<std::pair<edm::RefToBase<reco::Muon>, double>> &recAss = simToRec[sim];
    for (auto const &a : it->second) {
      recAss.push_back(std::make_pair(muons[a.idx], a.quality));
    }
  }
}