00001 00008 #include "FWCore/Framework/interface/EDProducer.h" 00009 #include "FWCore/ParameterSet/interface/InputTag.h" 00010 #include "DataFormats/Common/interface/Association.h" 00011 #include "DataFormats/HepMCCandidate/interface/GenParticle.h" 00012 00013 namespace edm { class ParameterSet; } 00014 00015 class MCTrackMatcher : public edm::EDProducer { 00016 public: 00018 MCTrackMatcher( const edm::ParameterSet & ); 00019 00020 private: 00021 void produce( edm::Event& evt, const edm::EventSetup& es ); 00022 std::string associator_; 00023 edm::InputTag tracks_, genParticles_, trackingParticles_; 00024 typedef edm::Association<reco::GenParticleCollection> GenParticleMatch; 00025 }; 00026 00027 #include "DataFormats/Common/interface/Handle.h" 00028 #include "FWCore/Framework/interface/ESHandle.h" 00029 #include "FWCore/Framework/interface/Event.h" 00030 #include "FWCore/Utilities/interface/EDMException.h" 00031 #include "FWCore/ParameterSet/interface/ParameterSet.h" 00032 #include "SimTracker/Records/interface/TrackAssociatorRecord.h" 00033 #include "SimTracker/TrackAssociation/interface/TrackAssociatorBase.h" 00034 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" 00035 using namespace edm; 00036 using namespace std; 00037 using namespace reco; 00038 00039 MCTrackMatcher::MCTrackMatcher(const ParameterSet & p) : 00040 associator_(p.getParameter<string>("associator")), 00041 tracks_(p.getParameter<InputTag>("tracks")), 00042 genParticles_( p.getParameter<InputTag>("genParticles")), 00043 trackingParticles_( p.getParameter<InputTag>("trackingParticles")) { 00044 produces<GenParticleMatch>(); 00045 } 00046 00047 void MCTrackMatcher::produce(Event& evt, const EventSetup& es) { 00048 ESHandle<TrackAssociatorBase> assoc; 00049 es.get<TrackAssociatorRecord>().get(associator_,assoc); 00050 const TrackAssociatorBase * associator = assoc.product(); 00051 Handle<View<Track> > tracks; 00052 evt.getByLabel(tracks_, tracks); 00053 Handle<TrackingParticleCollection> trackingParticles; 00054 evt.getByLabel(trackingParticles_,trackingParticles); 00055 Handle<vector<int> > barCodes; 00056 evt.getByLabel(genParticles_,barCodes ); 00057 Handle<GenParticleCollection> genParticles; 00058 evt.getByLabel(genParticles_, genParticles ); 00059 RecoToSimCollection associations = associator->associateRecoToSim ( tracks, trackingParticles, & evt ); 00060 auto_ptr<GenParticleMatch> match(new GenParticleMatch(GenParticleRefProd(genParticles))); 00061 GenParticleMatch::Filler filler(*match); 00062 size_t n = tracks->size(); 00063 vector<int> indices(n,-1); 00064 for (size_t i = 0; i < n; ++ i ) { 00065 RefToBase<Track> track(tracks, i); 00066 RecoToSimCollection::const_iterator f = associations.find(track); 00067 if ( f != associations.end() ) { 00068 TrackingParticleRef tp = f->val.front().first; 00069 const HepMC::GenParticle * particle = 0; 00070 TrackingParticle::genp_iterator j, b = tp->genParticle_begin(), e = tp->genParticle_end(); 00071 for( j = b; j != e; ++ j ) { 00072 const HepMC::GenParticle * p = j->get(); 00073 if (p->status() == 1) { 00074 particle = p; break; 00075 } 00076 } 00077 if( particle != 0 ) { 00078 int barCode = particle->barcode(); 00079 vector<int>::const_iterator 00080 b = barCodes->begin(), e = barCodes->end(), f = find( b, e, barCode ); 00081 if(f == e) throw edm::Exception(errors::InvalidReference) 00082 << "found matching particle with barcode" << *f 00083 << " which has not been found in " << genParticles_; 00084 indices[i] = *f; 00085 } 00086 } 00087 } 00088 filler.insert(tracks, indices.begin(), indices.end()); 00089 filler.fill(); 00090 evt.put(match); 00091 } 00092 00093 #include "FWCore/Framework/interface/MakerMacros.h" 00094 00095 DEFINE_FWK_MODULE( MCTrackMatcher ); 00096