Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-18 02:29:46

0001 #include "DataFormats/Candidate/interface/Candidate.h"
0002 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0003 #include "DataFormats/Common/interface/Association.h"
0004 #include "DataFormats/Common/interface/ValueMap.h"
0005 #include "DataFormats/Common/interface/View.h"
0006 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0007 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0008 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
0009 #include "DataFormats/Math/interface/deltaR.h"
0010 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0011 #include "DataFormats/PatCandidates/interface/Electron.h"
0012 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0013 #include "DataFormats/PatCandidates/interface/Photon.h"
0014 #include "DataFormats/VertexReco/interface/Vertex.h"
0015 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/Framework/interface/stream/EDProducer.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 
0023 #include <memory>
0024 
0025 // ------------------------------------------------------------------------------------------
0026 class PuppiPhoton : public edm::stream::EDProducer<> {
0027 public:
0028   explicit PuppiPhoton(const edm::ParameterSet &);
0029   ~PuppiPhoton() override;
0030 
0031   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0032   typedef math::XYZTLorentzVector LorentzVector;
0033   typedef std::vector<LorentzVector> LorentzVectorCollection;
0034   typedef edm::View<reco::Candidate> CandidateView;
0035   typedef std::vector<reco::PFCandidate> PFOutputCollection;
0036   typedef edm::View<reco::PFCandidate> PFView;
0037 
0038 private:
0039   void produce(edm::Event &, const edm::EventSetup &) override;
0040   bool matchPFCandidate(const reco::Candidate *iPF, const reco::Candidate *iPho);
0041   edm::EDGetTokenT<CandidateView> tokenPFCandidates_;
0042   edm::EDGetTokenT<CandidateView> tokenPuppiCandidates_;
0043   edm::EDGetTokenT<CandidateView> tokenPhotonCandidates_;
0044   edm::EDGetTokenT<edm::ValueMap<std::vector<reco::PFCandidateRef>>> reco2pf_;
0045   edm::EDGetTokenT<edm::ValueMap<float>> tokenWeights_;
0046   edm::EDGetTokenT<edm::ValueMap<bool>> tokenPhotonId_;
0047   double pt_;
0048   double eta_;
0049   bool usePFRef_;
0050   bool usePFphotons_;
0051   bool runOnMiniAOD_;
0052   bool usePhotonId_;
0053   std::vector<double> dRMatch_;
0054   std::vector<int32_t> pdgIds_;
0055   std::unique_ptr<PFOutputCollection> corrCandidates_;
0056   double weight_;
0057   bool useValueMap_;
0058 };
0059 
0060 // ------------------------------------------------------------------------------------------
0061 PuppiPhoton::PuppiPhoton(const edm::ParameterSet &iConfig) {
0062   tokenPFCandidates_ = consumes<CandidateView>(iConfig.getParameter<edm::InputTag>("candName"));
0063   tokenPuppiCandidates_ = consumes<CandidateView>(iConfig.getParameter<edm::InputTag>("puppiCandName"));
0064   usePFphotons_ = iConfig.getParameter<bool>("usePFphotons");
0065   if (!usePFphotons_)
0066     tokenPhotonCandidates_ = consumes<CandidateView>(iConfig.getParameter<edm::InputTag>("photonName"));
0067   usePhotonId_ = !(iConfig.getParameter<edm::InputTag>("photonId")).label().empty();
0068   if (usePhotonId_)
0069     tokenPhotonId_ = consumes<edm::ValueMap<bool>>(iConfig.getParameter<edm::InputTag>("photonId"));
0070   runOnMiniAOD_ = iConfig.getParameter<bool>("runOnMiniAOD");
0071   if (!runOnMiniAOD_)
0072     reco2pf_ =
0073         consumes<edm::ValueMap<std::vector<reco::PFCandidateRef>>>(iConfig.getParameter<edm::InputTag>("recoToPFMap"));
0074   useValueMap_ = iConfig.getParameter<bool>("useValueMap");
0075   if (useValueMap_)
0076     tokenWeights_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("weightsName"));
0077 
0078   pt_ = iConfig.getParameter<double>("pt");
0079   eta_ = iConfig.getParameter<double>("eta");
0080   dRMatch_ = iConfig.getParameter<std::vector<double>>("dRMatch");
0081   pdgIds_ = iConfig.getParameter<std::vector<int32_t>>("pdgids");
0082   usePFRef_ = iConfig.getParameter<bool>("useRefs");
0083   weight_ = iConfig.getParameter<double>("weight");
0084   produces<PFOutputCollection>();
0085   produces<edm::ValueMap<reco::CandidatePtr>>();
0086 }
0087 // ------------------------------------------------------------------------------------------
0088 PuppiPhoton::~PuppiPhoton() {}
0089 // ------------------------------------------------------------------------------------------
0090 void PuppiPhoton::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0091   int iC = -1;
0092   std::vector<const reco::Candidate *> phoCands;
0093   std::vector<uint16_t> phoIndx;
0094 
0095   edm::Handle<edm::ValueMap<std::vector<reco::PFCandidateRef>>> reco2pf;
0096   if (!runOnMiniAOD_)
0097     iEvent.getByToken(reco2pf_, reco2pf);
0098 
0099   // Get PFCandidate Collection
0100   edm::Handle<CandidateView> hPFProduct;
0101   iEvent.getByToken(tokenPFCandidates_, hPFProduct);
0102   const CandidateView *pfCol = hPFProduct.product();
0103 
0104   edm::Handle<CandidateView> hPuppiProduct;
0105   iEvent.getByToken(tokenPuppiCandidates_, hPuppiProduct);
0106   const CandidateView *pupCol = hPuppiProduct.product();
0107   if (usePFphotons_) {
0108     for (const auto &pho : *pfCol) {
0109       iC++;
0110       if (pho.pt() < pt_)
0111         continue;
0112       if (std::abs(pho.pdgId()) != 22)
0113         continue;
0114       if (fabs(pho.eta()) < eta_) {
0115         phoIndx.push_back(iC);
0116         phoCands.push_back(&pho);
0117       }
0118     }
0119   } else {
0120     edm::Handle<CandidateView> hPhoProduct;
0121     iEvent.getByToken(tokenPhotonCandidates_, hPhoProduct);
0122     const CandidateView *phoCol = hPhoProduct.product();
0123 
0124     edm::Handle<edm::ValueMap<bool>> photonId;
0125     if (usePhotonId_)
0126       iEvent.getByToken(tokenPhotonId_, photonId);
0127 
0128     for (CandidateView::const_iterator itPho = phoCol->begin(); itPho != phoCol->end(); itPho++) {
0129       iC++;
0130       bool passObject = false;
0131       if (itPho->isPhoton() && usePhotonId_)
0132         passObject = (*photonId)[phoCol->ptrAt(iC)];
0133       if (itPho->pt() < pt_)
0134         continue;
0135       if (!passObject && usePhotonId_)
0136         continue;
0137       if (runOnMiniAOD_) {
0138         const pat::Photon *pPho = dynamic_cast<const pat::Photon *>(&(*itPho));
0139         if (pPho != nullptr) {
0140           for (const edm::Ref<pat::PackedCandidateCollection> &ref : pPho->associatedPackedPFCandidates()) {
0141             if (fabs(ref->eta()) < eta_) {
0142               phoIndx.push_back(ref.key());
0143               phoCands.push_back(&(*(pfCol->ptrAt(ref.key()))));
0144             }
0145           }
0146           continue;
0147         }
0148         const pat::Electron *pElectron = dynamic_cast<const pat::Electron *>(&(*itPho));
0149         if (pElectron != nullptr) {
0150           for (const edm::Ref<pat::PackedCandidateCollection> &ref : pElectron->associatedPackedPFCandidates())
0151             if (fabs(ref->eta()) < eta_) {
0152               phoIndx.push_back(ref.key());
0153               phoCands.push_back(&(*(pfCol->ptrAt(ref.key()))));
0154             }
0155         }
0156       } else {
0157         for (const edm::Ref<std::vector<reco::PFCandidate>> &ref : (*reco2pf)[phoCol->ptrAt(iC)]) {
0158           if (fabs(ref->eta()) < eta_) {
0159             phoIndx.push_back(ref.key());
0160             phoCands.push_back(&(*(pfCol->ptrAt(ref.key()))));
0161           }
0162         }
0163       }
0164     }
0165   }
0166   //Get Weights
0167   edm::Handle<edm::ValueMap<float>> pupWeights;
0168   if (useValueMap_)
0169     iEvent.getByToken(tokenWeights_, pupWeights);
0170   std::unique_ptr<edm::ValueMap<LorentzVector>> p4PupOut(new edm::ValueMap<LorentzVector>());
0171   LorentzVectorCollection puppiP4s;
0172   std::vector<reco::CandidatePtr> values(hPFProduct->size());
0173   int iPF = 0;
0174   std::vector<float> lWeights;
0175   static const reco::PFCandidate dummySinceTranslateIsNotStatic;
0176   corrCandidates_ = std::make_unique<PFOutputCollection>();
0177   std::set<int> foundPhoIndex;
0178   for (CandidateView::const_iterator itPF = pupCol->begin(); itPF != pupCol->end(); itPF++) {
0179     auto id = dummySinceTranslateIsNotStatic.translatePdgIdToType(itPF->pdgId());
0180     const reco::PFCandidate *pPF = dynamic_cast<const reco::PFCandidate *>(&(*itPF));
0181     reco::PFCandidate pCand(pPF ? *pPF : reco::PFCandidate(itPF->charge(), itPF->p4(), id));
0182     LorentzVector pVec = itPF->p4();
0183     float pWeight = 1.;
0184     if (useValueMap_)
0185       pWeight = (*pupWeights)[pupCol->ptrAt(iPF)];
0186     if (!usePFRef_) {
0187       int iPho = -1;
0188       for (std::vector<const reco::Candidate *>::iterator itPho = phoCands.begin(); itPho != phoCands.end(); itPho++) {
0189         iPho++;
0190         if ((!matchPFCandidate(&(*itPF), *itPho)) || (foundPhoIndex.count(iPho) != 0))
0191           continue;
0192         pWeight = weight_;
0193         if (!useValueMap_ && itPF->pt() != 0)
0194           pWeight = pWeight * (phoCands[iPho]->pt() / itPF->pt());
0195         if (!useValueMap_ && itPF->pt() == 0)
0196           pVec.SetPxPyPzE(phoCands[iPho]->px() * pWeight,
0197                           phoCands[iPho]->py() * pWeight,
0198                           phoCands[iPho]->pz() * pWeight,
0199                           phoCands[iPho]->energy() * pWeight);
0200         foundPhoIndex.insert(iPho);
0201       }
0202     } else {
0203       int iPho = -1;
0204       for (std::vector<uint16_t>::const_iterator itPho = phoIndx.begin(); itPho != phoIndx.end(); itPho++) {
0205         iPho++;
0206         if (pupCol->refAt(iPF).key() != *itPho)
0207           continue;
0208         pWeight = weight_;
0209         if (!useValueMap_ && itPF->pt() != 0)
0210           pWeight = pWeight * (phoCands[iPho]->pt() / itPF->pt());
0211         if (!useValueMap_ && itPF->pt() == 0)
0212           pVec.SetPxPyPzE(phoCands[iPho]->px() * pWeight,
0213                           phoCands[iPho]->py() * pWeight,
0214                           phoCands[iPho]->pz() * pWeight,
0215                           phoCands[iPho]->energy() * pWeight);
0216         foundPhoIndex.insert(iPho);
0217       }
0218     }
0219     if (itPF->pt() != 0)
0220       pVec.SetPxPyPzE(itPF->px() * pWeight, itPF->py() * pWeight, itPF->pz() * pWeight, itPF->energy() * pWeight);
0221 
0222     lWeights.push_back(pWeight);
0223     pCand.setP4(pVec);
0224     puppiP4s.push_back(pVec);
0225     pCand.setSourceCandidatePtr(itPF->sourceCandidatePtr(0));
0226     corrCandidates_->push_back(pCand);
0227     iPF++;
0228   }
0229   //Add the missing pfcandidates
0230   for (unsigned int iPho = 0; iPho < phoCands.size(); iPho++) {
0231     if (foundPhoIndex.count(iPho) != 0)
0232       continue;
0233     auto id = dummySinceTranslateIsNotStatic.translatePdgIdToType(phoCands[iPho]->pdgId());
0234     reco::PFCandidate pCand(reco::PFCandidate(phoCands[iPho]->charge(), phoCands[iPho]->p4(), id));
0235     pCand.setSourceCandidatePtr(phoCands[iPho]->sourceCandidatePtr(0));
0236     LorentzVector pVec = phoCands[iPho]->p4();
0237     pVec.SetPxPyPzE(phoCands[iPho]->px() * weight_,
0238                     phoCands[iPho]->py() * weight_,
0239                     phoCands[iPho]->pz() * weight_,
0240                     phoCands[iPho]->energy() * weight_);
0241     pCand.setP4(pVec);
0242     lWeights.push_back(weight_);
0243     puppiP4s.push_back(pVec);
0244     corrCandidates_->push_back(pCand);
0245   }
0246   //Fill it into the event
0247   edm::OrphanHandle<reco::PFCandidateCollection> oh = iEvent.put(std::move(corrCandidates_));
0248   for (unsigned int ic = 0, nc = pupCol->size(); ic < nc; ++ic) {
0249     reco::CandidatePtr pkref(oh, ic);
0250     values[ic] = pkref;
0251   }
0252   std::unique_ptr<edm::ValueMap<reco::CandidatePtr>> pfMap_p(new edm::ValueMap<reco::CandidatePtr>());
0253   edm::ValueMap<reco::CandidatePtr>::Filler filler(*pfMap_p);
0254   filler.insert(hPFProduct, values.begin(), values.end());
0255   filler.fill();
0256   iEvent.put(std::move(pfMap_p));
0257 }
0258 // ------------------------------------------------------------------------------------------
0259 bool PuppiPhoton::matchPFCandidate(const reco::Candidate *iPF, const reco::Candidate *iPho) {
0260   if (iPF->pdgId() != iPho->pdgId())
0261     return false;
0262   double lDR = deltaR(iPF->eta(), iPF->phi(), iPho->eta(), iPho->phi());
0263   for (unsigned int i0 = 0; i0 < pdgIds_.size(); i0++) {
0264     if (std::abs(iPF->pdgId()) == pdgIds_[i0] && lDR < dRMatch_[i0])
0265       return true;
0266   }
0267   return false;
0268 }
0269 // ------------------------------------------------------------------------------------------
0270 void PuppiPhoton::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0271   //The following says we do not know what parameters are allowed so do no validation
0272   // Please change this to state exactly what you do use, even if it is no parameters
0273   edm::ParameterSetDescription desc;
0274   desc.setUnknown();
0275   descriptions.addDefault(desc);
0276 }
0277 //define this as a plug-in
0278 DEFINE_FWK_MODULE(PuppiPhoton);