NanoAOD/plugins/BJetEnergyRegressionMVA.cc

0001 //
0002 // Original Author:  Andre Rizzi
0003 //         Created:  Mon, 07 Sep 2017 09:18:03 GMT
0004 //
0005 //
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Utilities/interface/StreamID.h"
0015
0016 #include "DataFormats/PatCandidates/interface/Muon.h"
0017 #include "DataFormats/PatCandidates/interface/Electron.h"
0018 #include "DataFormats/PatCandidates/interface/Jet.h"
0019
0020 #include "DataFormats/VertexReco/interface/Vertex.h"
0021 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0022
0023 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
0024 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0025 #include "RecoVertex/VertexPrimitives/interface/VertexState.h"
0026
0027 #include "PhysicsTools/PatAlgos/interface/BaseMVAValueMapProducer.h"
0028 #include <vector>
0029
0030 class BJetEnergyRegressionMVA : public BaseMVAValueMapProducer<pat::Jet> {
0031 public:
0032   explicit BJetEnergyRegressionMVA(const edm::ParameterSet& iConfig, const BaseMVACache* cache)
0033       : BaseMVAValueMapProducer<pat::Jet>(iConfig, cache),
0034         pvsrc_(consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvsrc"))),
0035         svsrc_(consumes<edm::View<reco::VertexCompositePtrCandidate>>(iConfig.getParameter<edm::InputTag>("svsrc"))),
0036         rhosrc_(consumes<double>(iConfig.getParameter<edm::InputTag>("rhosrc"))) {}
0037
0038   void readAdditionalCollections(edm::Event& iEvent, const edm::EventSetup&) override {
0039     iEvent.getByToken(pvsrc_, pvs_);
0040     iEvent.getByToken(svsrc_, svs_);
0041     iEvent.getByToken(rhosrc_, rho_);
0042   }
0043
0044   void fillAdditionalVariables(const pat::Jet& j) override {
0045     this->setValue("nPVs", pvs_->size());
0046     this->setValue("rho", *(rho_.product()));
0047
0048     float cone_boundaries[] = {0.05, 0.1, 0.2, 0.3, 0.4};
0049     size_t ncone_boundaries = sizeof(cone_boundaries) / sizeof(float);
0050     std::vector<float> emFractionEnergyRings(ncone_boundaries + 1);
0051     std::vector<float> chFractionEnergyRings(ncone_boundaries + 1);
0052     std::vector<float> neFractionEnergyRings(ncone_boundaries + 1);
0053     std::vector<float> muFractionEnergyRings(ncone_boundaries + 1);
0054     float jetRawEnergy = j.p4().E() * j.jecFactor("Uncorrected");
0055     int numDaughtersPt03 = 0;
0056     for (unsigned int ijcone = 0; ijcone < ncone_boundaries; ijcone++) {
0057       emFractionEnergyRings[ijcone] = 0;
0058       muFractionEnergyRings[ijcone] = 0;
0059       chFractionEnergyRings[ijcone] = 0;
0060       neFractionEnergyRings[ijcone] = 0;
0061     }
0062     for (const auto& d : j.daughterPtrVector()) {
0063       float candDr = Geom::deltaR(d->p4(), j.p4());
0064       size_t icone =
0065           std::lower_bound(&cone_boundaries[0], &cone_boundaries[ncone_boundaries], candDr) - &cone_boundaries[0];
0066       float candEnergy = d->energy() / jetRawEnergy;
0067       int pdgid = abs(d->pdgId());
0068       if (pdgid == 22 || pdgid == 11) {
0069         emFractionEnergyRings[icone] += candEnergy;
0070       } else if (pdgid == 13) {
0071         muFractionEnergyRings[icone] += candEnergy;
0072       } else if (d->charge() != 0) {
0073         chFractionEnergyRings[icone] += candEnergy;
0074       } else {
0075         neFractionEnergyRings[icone] += candEnergy;
0076       }
0077       if (d->pt() > 0.3)
0078         numDaughtersPt03 += 1;
0079     }  // end of jet daughters loop
0080
0081     this->setValue("Jet_energyRing_dR0_em_Jet_rawEnergy", emFractionEnergyRings[0]);
0082     this->setValue("Jet_energyRing_dR1_em_Jet_rawEnergy", emFractionEnergyRings[1]);
0083     this->setValue("Jet_energyRing_dR2_em_Jet_rawEnergy", emFractionEnergyRings[2]);
0084     this->setValue("Jet_energyRing_dR3_em_Jet_rawEnergy", emFractionEnergyRings[3]);
0085     this->setValue("Jet_energyRing_dR4_em_Jet_rawEnergy", emFractionEnergyRings[4]);
0086     //  this->setValue("Jet_energyRing_dR5_em_Jet_rawEnergy", emFractionEnergyRings[5]);
0087
0088     this->setValue("Jet_energyRing_dR0_ch_Jet_rawEnergy", chFractionEnergyRings[0]);
0089     this->setValue("Jet_energyRing_dR1_ch_Jet_rawEnergy", chFractionEnergyRings[1]);
0090     this->setValue("Jet_energyRing_dR2_ch_Jet_rawEnergy", chFractionEnergyRings[2]);
0091     this->setValue("Jet_energyRing_dR3_ch_Jet_rawEnergy", chFractionEnergyRings[3]);
0092     this->setValue("Jet_energyRing_dR4_ch_Jet_rawEnergy", chFractionEnergyRings[4]);
0093     //  this->setValue("Jet_energyRing_dR5_ch_Jet_rawEnergy", chFractionEnergyRings[5]);
0094
0095     this->setValue("Jet_energyRing_dR0_mu_Jet_rawEnergy", muFractionEnergyRings[0]);
0096     this->setValue("Jet_energyRing_dR1_mu_Jet_rawEnergy", muFractionEnergyRings[1]);
0097     this->setValue("Jet_energyRing_dR2_mu_Jet_rawEnergy", muFractionEnergyRings[2]);
0098     this->setValue("Jet_energyRing_dR3_mu_Jet_rawEnergy", muFractionEnergyRings[3]);
0099     this->setValue("Jet_energyRing_dR4_mu_Jet_rawEnergy", muFractionEnergyRings[4]);
0100     //  this->setValue("Jet_energyRing_dR5_mu_Jet_rawEnergy", muFractionEnergyRings[5]);
0101
0102     this->setValue("Jet_energyRing_dR0_neut_Jet_rawEnergy", neFractionEnergyRings[0]);
0103     this->setValue("Jet_energyRing_dR1_neut_Jet_rawEnergy", neFractionEnergyRings[1]);
0104     this->setValue("Jet_energyRing_dR2_neut_Jet_rawEnergy", neFractionEnergyRings[2]);
0105     this->setValue("Jet_energyRing_dR3_neut_Jet_rawEnergy", neFractionEnergyRings[3]);
0106     this->setValue("Jet_energyRing_dR4_neut_Jet_rawEnergy", neFractionEnergyRings[4]);
0107     //    this->setValue("Jet_energyRing_dR5_neut_Jet_rawEnergy", neFractionEnergyRings[5]);
0108
0109     this->setValue("Jet_numDaughters_pt03", numDaughtersPt03);
0110   }
0111
0112   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0113     edm::ParameterSetDescription desc = BaseMVAValueMapProducer<pat::Jet>::getDescription();
0114     desc.add<edm::InputTag>("pvsrc")->setComment("primary vertices input collection");
0115     desc.add<edm::InputTag>("svsrc")->setComment("secondary vertices input collection");
0116     desc.add<edm::InputTag>("rhosrc")->setComment("rho  input collection");
0117     descriptions.add("BJetEnergyRegressionMVA", desc);
0118   }
0119
0120 private:
0121   const edm::EDGetTokenT<std::vector<reco::Vertex>> pvsrc_;
0122   edm::Handle<std::vector<reco::Vertex>> pvs_;
0123   const edm::EDGetTokenT<edm::View<reco::VertexCompositePtrCandidate>> svsrc_;
0124   edm::Handle<edm::View<reco::VertexCompositePtrCandidate>> svs_;
0125   edm::EDGetTokenT<double> rhosrc_;
0126   edm::Handle<double> rho_;
0127 };
0128
0129 //define this as a plug-in
0130 DEFINE_FWK_MODULE(BJetEnergyRegressionMVA);