Go to the documentation of this file.00001
00008 #include "FWCore/Framework/interface/EDProducer.h"
00009 #include "FWCore/Utilities/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