MuonAssociatorByHits.cc

Go to the documentation of this file.

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

using namespace reco;
using namespace std;
using namespace muonAssociatorByHitsDiagnostics;

namespace muonAssociatorByHitsDiagnostics {
  using TrackHitsCollection = MuonAssociatorByHitsHelper::TrackHitsCollection;

  class InputDumper {
  public:
    InputDumper(const edm::ParameterSet& conf) :
      simtracksTag(conf.getParameter<edm::InputTag>("simtracksTag")),
      simtracksXFTag(conf.getParameter<edm::InputTag>("simtracksXFTag")),
      crossingframe(conf.getParameter<bool>("crossingframe")) {}

    InputDumper(const edm::ParameterSet& conf, edm::ConsumesCollector && iC) :
      InputDumper(conf)
    {
      if (crossingframe) {
        iC.consumes<CrossingFrame<SimTrack> >(simtracksXFTag);
        iC.consumes<CrossingFrame<SimVertex> >(simtracksXFTag);
      }
      else{
        iC.consumes<edm::SimTrackContainer>(simtracksTag);
        iC.consumes<edm::SimVertexContainer>(simtracksTag);
      }
    }

    void dump(const TrackHitsCollection&, const TrackingParticleCollection&, const edm::Event&) const;

  private:
    edm::InputTag const simtracksTag;
    edm::InputTag const simtracksXFTag;
    bool const crossingframe;
  };

  void
  InputDumper::dump(const TrackHitsCollection& tC, const TrackingParticleCollection& tPC, const edm::Event& event) const {
    // reco::Track collection
    edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"reco::Track collection --- size = "<<tC.size();
 
 
    // TrackingParticle collection
    edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"TrackingParticle collection --- size = "<<tPC.size();
    int j = 0;
    for(TrackingParticleCollection::const_iterator ITER=tPC.begin(); ITER!=tPC.end(); ITER++, j++) {
      edm::LogVerbatim("MuonAssociatorByHits")
        <<"TrackingParticle "<<j<<", q = "<<ITER->charge()<<", p = "<<ITER->p()
        <<", pT = "<<ITER->pt()<<", eta = "<<ITER->eta()<<", phi = "<<ITER->phi();
       
      edm::LogVerbatim("MuonAssociatorByHits")
        <<"\t pdg code = "<<ITER->pdgId()<<", made of "<<ITER->numberOfHits()<<" PSimHit"
        <<" (in "<<ITER->numberOfTrackerLayers()<<" layers)"
        <<" from "<<ITER->g4Tracks().size()<<" SimTrack:";
      for (TrackingParticle::g4t_iterator g4T=ITER->g4Track_begin(); g4T!=ITER->g4Track_end(); g4T++) {
        edm::LogVerbatim("MuonAssociatorByHits")
          <<"\t\t Id:"<<g4T->trackId()<<"/Evt:("<<g4T->eventId().event()<<","<<g4T->eventId().bunchCrossing()<<")";
      }    
    }
 
    // SimTrack collection
    edm::Handle<CrossingFrame<SimTrack> > cf_simtracks;
    edm::Handle<edm::SimTrackContainer> simTrackCollection;
 
    // SimVertex collection
    edm::Handle<CrossingFrame<SimVertex> > cf_simvertices;
    edm::Handle<edm::SimVertexContainer> simVertexCollection;
         
    if (crossingframe) {
      event.getByLabel(simtracksXFTag,cf_simtracks);
      auto_ptr<MixCollection<SimTrack> > SimTk( new MixCollection<SimTrack>(cf_simtracks.product()) );
      edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"CrossingFrame<SimTrack> collection with InputTag = "<<simtracksXFTag
                                              <<" has size = "<<SimTk->size();
      int k = 0;
      for (MixCollection<SimTrack>::MixItr ITER=SimTk->begin(); ITER!=SimTk->end(); ITER++, k++) {
        edm::LogVerbatim("MuonAssociatorByHits")
          <<"SimTrack "<<k
          <<" - Id:"<<ITER->trackId()<<"/Evt:("<<ITER->eventId().event()<<","<<ITER->eventId().bunchCrossing()<<")"
          <<" pdgId = "<<ITER->type()<<", q = "<<ITER->charge()<<", p = "<<ITER->momentum().P()
          <<", pT = "<<ITER->momentum().Pt()<<", eta = "<<ITER->momentum().Eta()<<", phi = "<<ITER->momentum().Phi()
          <<"\n * "<<*ITER <<endl;
      }
      event.getByLabel(simtracksXFTag,cf_simvertices);
      auto_ptr<MixCollection<SimVertex> > SimVtx( new MixCollection<SimVertex>(cf_simvertices.product()) );
      edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"CrossingFrame<SimVertex> collection with InputTag = "<<simtracksXFTag
                                              <<" has size = "<<SimVtx->size();
      int kv = 0;
      for (MixCollection<SimVertex>::MixItr VITER=SimVtx->begin(); VITER!=SimVtx->end(); VITER++, kv++){
        edm::LogVerbatim("MuonAssociatorByHits")
          <<"SimVertex "<<kv
          << " : "<< *VITER <<endl;
      }
    }
    else {
      event.getByLabel(simtracksTag,simTrackCollection);
      const edm::SimTrackContainer simTC = *(simTrackCollection.product());
      edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"SimTrack collection with InputTag = "<<simtracksTag
                                              <<" has size = "<<simTC.size()<<endl;
      int k = 0;
      for(edm::SimTrackContainer::const_iterator ITER=simTC.begin(); ITER!=simTC.end(); ITER++, k++){
        edm::LogVerbatim("MuonAssociatorByHits")
          <<"SimTrack "<<k
          <<" - Id:"<<ITER->trackId()<<"/Evt:("<<ITER->eventId().event()<<","<<ITER->eventId().bunchCrossing()<<")"
          <<" pdgId = "<<ITER->type()<<", q = "<<ITER->charge()<<", p = "<<ITER->momentum().P()
          <<", pT = "<<ITER->momentum().Pt()<<", eta = "<<ITER->momentum().Eta()<<", phi = "<<ITER->momentum().Phi()
          <<"\n * "<<*ITER <<endl;
      }
      event.getByLabel(simtracksTag,simVertexCollection);
      const edm::SimVertexContainer simVC = *(simVertexCollection.product());
      edm::LogVerbatim("MuonAssociatorByHits")<<"\n"<<"SimVertex collection with InputTag = "<<"g4SimHits"
                                              <<" has size = "<<simVC.size()<<endl;
      int kv = 0;
      for (edm::SimVertexContainer::const_iterator VITER=simVC.begin(); VITER!=simVC.end(); VITER++, kv++){
        edm::LogVerbatim("MuonAssociatorByHits")
          <<"SimVertex "<<kv
          << " : "<< *VITER <<endl;
      }
    }
  }

}


MuonAssociatorByHits::MuonAssociatorByHits (const edm::ParameterSet& conf, edm::ConsumesCollector && iC) :  
  helper_(conf),
  conf_(conf),
  trackerHitAssociatorConfig_(conf, std::move(iC))
{
  //hack for consumes
  RPCHitAssociator rpctruth(conf,std::move(iC));
  DTHitAssociator dttruth(conf,std::move(iC));
  CSCHitAssociator muonTruth(conf,std::move(iC));
  if( conf.getUntrackedParameter<bool>("dumpInputCollections") ) {
    diagnostics_.reset( new InputDumper(conf, std::move(iC)) );
  }
}



MuonAssociatorByHits::~MuonAssociatorByHits()
{
}

RecoToSimCollection  
MuonAssociatorByHits::associateRecoToSim( const edm::RefToBaseVector<reco::Track>& tC,
                      const edm::RefVector<TrackingParticleCollection>& TPCollectionH,
                          const edm::Event * e, const edm::EventSetup * setup) const{
  RecoToSimCollection  outputCollection(&e->productGetter());

  MuonAssociatorByHitsHelper::TrackHitsCollection tH;
  for (auto it = tC.begin(), ed = tC.end(); it != ed; ++it) {
    tH.push_back(std::make_pair((*it)->recHitsBegin(), (*it)->recHitsEnd()));
  }

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


  // Tracker hit association  
  TrackerHitAssociator trackertruth(*e, trackerHitAssociatorConfig_);
  // CSC hit association
  CSCHitAssociator csctruth(*e,*setup,conf_);
  // DT hit association
  bool printRtS(true);
  DTHitAssociator dttruth(*e,*setup,conf_,printRtS);  
  // RPC hit association
  RPCHitAssociator rpctruth(*e,*setup,conf_);
   
  MuonAssociatorByHitsHelper::Resources resources = {tTopo, &trackertruth, &csctruth, &dttruth, &rpctruth};

  if(diagnostics_) {
    resources.diagnostics_ = [this, e](const TrackHitsCollection& hC, const TrackingParticleCollection& pC) {
      diagnostics_->dump(hC,pC, *e);
    };
  }

  auto bareAssoc = helper_.associateRecoToSimIndices(tH, TPCollectionH, resources);
  for (auto it = bareAssoc.begin(), ed = bareAssoc.end(); it != ed; ++it) {
    for (auto itma = it->second.begin(), edma = it->second.end(); itma != edma; ++itma) {
        outputCollection.insert(tC[it->first], std::make_pair(TPCollectionH[itma->idx], itma->quality));
    }
  }

  outputCollection.post_insert(); // perhaps not even necessary
  return outputCollection;
}

SimToRecoCollection  
MuonAssociatorByHits::associateSimToReco( const edm::RefToBaseVector<reco::Track>& tC, 
                      const edm::RefVector<TrackingParticleCollection>& TPCollectionH,
                      const edm::Event * e, const edm::EventSetup * setup) const{

  SimToRecoCollection  outputCollection(&e->productGetter());
  MuonAssociatorByHitsHelper::TrackHitsCollection tH;
  for (auto it = tC.begin(), ed = tC.end(); it != ed; ++it) {
    tH.push_back(std::make_pair((*it)->recHitsBegin(), (*it)->recHitsEnd()));
  }

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

  // Tracker hit association  
  TrackerHitAssociator trackertruth(*e, trackerHitAssociatorConfig_);
  // CSC hit association
  CSCHitAssociator csctruth(*e,*setup,conf_);
  // DT hit association
  bool printRtS = false;
  DTHitAssociator dttruth(*e,*setup,conf_,printRtS);  
  // RPC hit association
  RPCHitAssociator rpctruth(*e,*setup,conf_);
   
  MuonAssociatorByHitsHelper::Resources resources = {tTopo, &trackertruth, &csctruth, &dttruth, &rpctruth};
  
  auto bareAssoc = helper_.associateSimToRecoIndices(tH, TPCollectionH, resources);
  for (auto it = bareAssoc.begin(), ed = bareAssoc.end(); it != ed; ++it) {
    for (auto itma = it->second.begin(), edma = it->second.end(); itma != edma; ++itma) {
        outputCollection.insert(TPCollectionH[ it->first],
                                std::make_pair(tC[itma->idx], itma->quality));
    }
  }

  outputCollection.post_insert(); // perhaps not even necessary
  return outputCollection;
}