/ElectroWeakAnalysis/ZMuMu/plugins/DimuonMCMatcher.cc

 1#include "FWCore/Framework/interface/EDProducer.h"
 2#include "FWCore/Utilities/interface/InputTag.h"
 3#include <iostream>
 4
 5class DimuonMCMatcher : public edm::EDProducer {
 6public:
 7  DimuonMCMatcher(const edm::ParameterSet & cfg);
 8  virtual void produce(edm::Event&, const edm::EventSetup&) override;
 9private:
10  edm::InputTag src_;
11};
12
13#include "DataFormats/Candidate/interface/Candidate.h"
14#include "DataFormats/Candidate/interface/CandidateFwd.h"
15#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
16#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
17#include "FWCore/Framework/interface/Event.h"
18#include "DataFormats/Common/interface/Handle.h"
19#include "FWCore/ParameterSet/interface/ParameterSet.h"
20#include "DataFormats/PatCandidates/interface/Muon.h"
21#include "DataFormats/PatCandidates/interface/GenericParticle.h"
22using namespace std;
23using namespace reco;
24using namespace edm;
25
26DimuonMCMatcher::DimuonMCMatcher(const edm::ParameterSet & cfg) : 
27  src_(cfg.getParameter<InputTag>("src")) {
28  produces<vector<GenParticleRef> >();
29}
30
31void DimuonMCMatcher::produce(edm::Event& evt, const edm::EventSetup&) {
32  Handle<CandidateView> src;
33  evt.getByLabel(src_, src);
34  auto_ptr<vector<GenParticleRef> > matched(new vector<GenParticleRef>);
35  matched->reserve(src->size());
36  int j=0;
37  for(CandidateView::const_iterator i = src->begin(); i != src->end(); ++i) {
38    j++;
39    const Candidate * dau1 = i->daughter(0);
40    const Candidate * dau2 = i->daughter(1);
41    if(dau1 == 0|| dau2 == 0) 
42      throw Exception(errors::InvalidReference) <<
43	"one of the two daughter does not exist\n";
44    const Candidate * c1 = dau1->masterClone().get();
45    GenParticleRef mc1;
46    const pat::Muon * mu1 = dynamic_cast<const pat::Muon*>(c1);
47    if(mu1 != 0) {
48      mc1 = mu1->genParticleRef();
49      //     if (mc1.isNonnull()) cout << "DimuonMCMatcher> genParticleRef1 " << mc1->pdgId() << endl;
50    } else {
51      const pat::GenericParticle * gp1 = dynamic_cast<const pat::GenericParticle*>(c1);
52      if(gp1 == 0) 
53	throw Exception(errors::InvalidReference) <<
54	  "first of two daughter is neither a pat::Muon not pat::GenericParticle\n";
55      mc1 = gp1->genParticleRef();
56    }
57    const Candidate * c2 = dau2->masterClone().get();
58    GenParticleRef mc2;
59    const pat::Muon * mu2 = dynamic_cast<const pat::Muon*>(c2);
60    if(mu2 != 0) {
61      mc2 = mu2->genParticleRef();
62      //      if (mc2.isNonnull()) cout << "DimuonMCMatcher> genParticleRef2 " << mc2->pdgId() << endl;
63    } else {
64      const pat::GenericParticle * gp2 = dynamic_cast<const pat::GenericParticle*>(c2);
65      if(gp2 == 0) 
66	throw Exception(errors::InvalidReference) <<
67	  "first of two daughter is neither a pat::Muon not pat::GenericParticle\n";
68      mc2 = gp2->genParticleRef();
69    }
70    GenParticleRef dimuonMatch;
71    //    cout << "DimuonMatcher>  mc1 " << mc1.isNonnull() << "   mc2 " << mc2.isNonnull() << endl;
72    if(mc1.isNonnull() && mc2.isNonnull()) {
73      int k=0;
74      do {
75	k++;
76	mc1 = mc1->numberOfMothers() > 0 ? mc1->motherRef() : GenParticleRef();
77	mc2 = mc2->numberOfMothers() > 0 ? mc2->motherRef() : GenParticleRef();
78	//	cout << "DimuonMCMatcher> do loop: " << k << "  id1 " << mc1->pdgId() << "    id2 " << mc2->pdgId() << endl;
79      } while(mc1 != mc2 && mc1.isNonnull() && mc2.isNonnull());
80      if(mc1.isNonnull() && mc2.isNonnull() && mc1->pdgId()==23) {
81	dimuonMatch = mc1;
82      }
83    }
84    //    cout << "DimuonMatcher> dimuonMatch " << dimuonMatch.isNonnull() << endl;
85    matched->push_back(dimuonMatch);
86  }
87
88  evt.put(matched);
89}
90
91
92#include "FWCore/Framework/interface/MakerMacros.h"
93
94DEFINE_FWK_MODULE(DimuonMCMatcher);