JetProducers/plugins/QGTagger.cc

0001 #include <memory>
0002
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "DataFormats/VertexReco/interface/Vertex.h"
0010 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0011 #include "DataFormats/JetReco/interface/PFJet.h"
0012 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0013 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0014 #include "DataFormats/Common/interface/ValueMap.h"
0015 #include "DataFormats/Math/interface/deltaPhi.h"
0016
0017 #include "RecoJets/JetProducers/interface/QGTagger.h"
0018 #include "RecoJets/JetAlgorithms/interface/QGLikelihoodCalculator.h"
0019
0020 /**
0021  * EDProducer class to produced the qgLikelihood values and related variables
0022  * If the input jets are uncorrected, the jecService should be provided, so jet are corrected on the fly before the algorithm is applied
0023  * Authors: andrea.carlo.marini@cern.ch, tom.cornelis@cern.ch, cms-qg-workinggroup@cern.ch
0024  */
0025 QGTagger::QGTagger(const edm::ParameterSet& iConfig)
0026     : jetsToken(consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("srcJets"))),
0027       jetCorrectorToken(consumes<reco::JetCorrector>(iConfig.getParameter<edm::InputTag>("jec"))),
0028       vertexToken(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("srcVertexCollection"))),
0029       rhoToken(consumes<double>(iConfig.getParameter<edm::InputTag>("srcRho"))),
0030       computeLikelihood(iConfig.getParameter<bool>("computeLikelihood")),
0031       paramsToken(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("jetsLabel")))),
0032       useQC(iConfig.getParameter<bool>("useQualityCuts")),
0033       useJetCorr(!iConfig.getParameter<edm::InputTag>("jec").label().empty()),
0034       produceSyst(!iConfig.getParameter<std::string>("systematicsLabel").empty()),
0035       applyConstituentWeight(false) {
0036   produces<edm::ValueMap<float>>("axis2");
0037   produces<edm::ValueMap<int>>("mult");
0038   produces<edm::ValueMap<float>>("ptD");
0039   if (computeLikelihood) {
0040     produces<edm::ValueMap<float>>("qgLikelihood");
0041   }
0042
0043   edm::InputTag srcConstituentWeights = iConfig.getParameter<edm::InputTag>("srcConstituentWeights");
0044   if (!srcConstituentWeights.label().empty()) {
0045     constituentWeightsToken = consumes<edm::ValueMap<float>>(srcConstituentWeights);
0046     applyConstituentWeight = true;
0047   }
0048
0049   if (computeLikelihood && produceSyst) {
0050     systToken = esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("systematicsLabel")));
0051     produces<edm::ValueMap<float>>("qgLikelihoodSmearedQuark");
0052     produces<edm::ValueMap<float>>("qgLikelihoodSmearedGluon");
0053     produces<edm::ValueMap<float>>("qgLikelihoodSmearedAll");
0054   }
0055 }
0056
0057 /// Produce qgLikelihood using {mult, ptD, -log(axis2)}
0058 void QGTagger::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0059   std::vector<float> qgProduct;
0060   std::vector<float> axis2Product;
0061   std::vector<int> multProduct;
0062   std::vector<float> ptDProduct;
0063   std::vector<float> smearedQuarkProduct;
0064   std::vector<float> smearedGluonProduct;
0065   std::vector<float> smearedAllProduct;
0066
0067   edm::Handle<edm::View<reco::Jet>> jets;
0068   iEvent.getByToken(jetsToken, jets);
0069   edm::Handle<reco::JetCorrector> jetCorr;
0070   if (useJetCorr)
0071     iEvent.getByToken(jetCorrectorToken, jetCorr);
0072   edm::Handle<reco::VertexCollection> vertexCollection;
0073   iEvent.getByToken(vertexToken, vertexCollection);
0074   edm::Handle<double> rho;
0075   iEvent.getByToken(rhoToken, rho);
0076
0077   edm::ValueMap<float> constituentWeights;
0078   if (applyConstituentWeight) {
0079     constituentWeights = iEvent.get(constituentWeightsToken);
0080   }
0081
0082   const QGLikelihoodObject* QGLParamsColl = nullptr;
0083   if (computeLikelihood) {
0084     QGLParamsColl = &iSetup.getData(paramsToken);
0085   }
0086
0087   const QGLikelihoodSystematicsObject* QGLSystColl = nullptr;
0088   if (produceSyst) {
0089     QGLSystColl = &iSetup.getData(systToken);
0090   }
0091
0092   bool weStillNeedToCheckJetCandidates = true;
0093   bool weAreUsingPackedCandidates = false;
0094   for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
0095     if (weStillNeedToCheckJetCandidates) {
0096       weAreUsingPackedCandidates = isPackedCandidate(&*jet);
0097       weStillNeedToCheckJetCandidates = false;
0098     }
0099     float pt = (useJetCorr ? jet->pt() * jetCorr->correction(*jet) : jet->pt());
0100
0101     float ptD, axis2;
0102     int mult;
0103     std::tie(mult, ptD, axis2) = calcVariables(&*jet, vertexCollection, constituentWeights, weAreUsingPackedCandidates);
0104
0105     float qgValue;
0106     if (mult > 2 && computeLikelihood)
0107       qgValue =
0108           qgLikelihood.computeQGLikelihood(*QGLParamsColl, pt, jet->eta(), *rho, {(float)mult, ptD, -std::log(axis2)});
0109     else
0110       qgValue = -1;
0111
0112     qgProduct.push_back(qgValue);
0113     if (computeLikelihood && produceSyst) {
0114       smearedQuarkProduct.push_back(qgLikelihood.systematicSmearing(*QGLSystColl, pt, jet->eta(), *rho, qgValue, 0));
0115       smearedGluonProduct.push_back(qgLikelihood.systematicSmearing(*QGLSystColl, pt, jet->eta(), *rho, qgValue, 1));
0116       smearedAllProduct.push_back(qgLikelihood.systematicSmearing(*QGLSystColl, pt, jet->eta(), *rho, qgValue, 2));
0117     }
0118     axis2Product.push_back(axis2);
0119     multProduct.push_back(mult);
0120     ptDProduct.push_back(ptD);
0121   }
0122
0123   putInEvent("axis2", jets, axis2Product, iEvent);
0124   putInEvent("mult", jets, multProduct, iEvent);
0125   putInEvent("ptD", jets, ptDProduct, iEvent);
0126   if (computeLikelihood) {
0127     putInEvent("qgLikelihood", jets, qgProduct, iEvent);
0128     if (produceSyst) {
0129       putInEvent("qgLikelihoodSmearedQuark", jets, smearedQuarkProduct, iEvent);
0130       putInEvent("qgLikelihoodSmearedGluon", jets, smearedGluonProduct, iEvent);
0131       putInEvent("qgLikelihoodSmearedAll", jets, smearedAllProduct, iEvent);
0132     }
0133   }
0134 }
0135
0136 /// Function to put product into event
0137 template <typename T>
0138 void QGTagger::putInEvent(const std::string& name,
0139                           const edm::Handle<edm::View<reco::Jet>>& jets,
0140                           const std::vector<T>& product,
0141                           edm::Event& iEvent) const {
0142   auto out = std::make_unique<edm::ValueMap<T>>();
0143   typename edm::ValueMap<T>::Filler filler(*out);
0144   filler.insert(jets, product.begin(), product.end());
0145   filler.fill();
0146   iEvent.put(std::move(out), name);
0147 }
0148
0149 /// Function to tell us if we are using packedCandidates, only test for first candidate
0150 bool QGTagger::isPackedCandidate(const reco::Jet* jet) const {
0151   for (auto candidate : jet->getJetConstituentsQuick()) {
0152     if (typeid(pat::PackedCandidate) == typeid(*candidate))
0153       return true;
0154     else if (typeid(reco::PFCandidate) == typeid(*candidate))
0155       return false;
0156     else
0157       throw cms::Exception("WrongJetCollection", "Jet constituents are not particle flow candidates");
0158   }
0159   return false;
0160 }
0161
0162 /// Calculation of axis2, mult and ptD
0163 std::tuple<int, float, float> QGTagger::calcVariables(const reco::Jet* jet,
0164                                                       edm::Handle<reco::VertexCollection>& vC,
0165                                                       edm::ValueMap<float>& constituentWeights,
0166                                                       bool weAreUsingPackedCandidates) const {
0167   float sum_weight = 0., sum_deta = 0., sum_dphi = 0., sum_deta2 = 0., sum_dphi2 = 0., sum_detadphi = 0., sum_pt = 0.;
0168
0169   float multWeighted = 0;
0170
0171   //Loop over the jet constituents
0172   for (const auto& daughter : jet->getJetConstituents()) {
0173     const reco::Candidate* cand = daughter.get();
0174
0175     float constWeight = 1.0;
0176     if (applyConstituentWeight) {
0177       constWeight = constituentWeights[daughter];
0178     }
0179
0180     if (weAreUsingPackedCandidates) {  //packed candidate situation
0181       auto part = static_cast<const pat::PackedCandidate*>(cand);
0182
0183       if (part->charge()) {
0184         if (!(part->fromPV() > 1 && part->trackHighPurity()))
0185           continue;
0186         if (useQC) {
0187           if ((part->dz() * part->dz()) / (part->dzError() * part->dzError()) > 25.)
0188             continue;
0189           if ((part->dxy() * part->dxy()) / (part->dxyError() * part->dxyError()) < 25.)
0190             multWeighted += constWeight;
0191         } else
0192           multWeighted += constWeight;
0193       } else {
0194         if ((constWeight * part->pt()) < 1.0)
0195           continue;
0196         multWeighted += constWeight;
0197       }
0198     } else {
0199       auto part = static_cast<const reco::PFCandidate*>(cand);
0200
0201       reco::TrackRef itrk = part->trackRef();
0202       if (itrk.isNonnull()) {  //Track exists --> charged particle
0203         auto vtxLead = vC->begin();
0204         auto vtxClose = vC->begin();  //Search for closest vertex to track
0205         for (auto vtx = vC->begin(); vtx != vC->end(); ++vtx) {
0206           if (fabs(itrk->dz(vtx->position())) < fabs(itrk->dz(vtxClose->position())))
0207             vtxClose = vtx;
0208         }
0209         if (!(vtxClose == vtxLead && itrk->quality(reco::TrackBase::qualityByName("highPurity"))))
0210           continue;
0211
0212         if (useQC) {  //If useQC, require dz and d0 cuts
0213           float dz = itrk->dz(vtxClose->position());
0214           float d0 = itrk->dxy(vtxClose->position());
0215           float dz_sigma_square = pow(itrk->dzError(), 2) + pow(vtxClose->zError(), 2);
0216           float d0_sigma_square = pow(itrk->d0Error(), 2) + pow(vtxClose->xError(), 2) + pow(vtxClose->yError(), 2);
0217           if (dz * dz / dz_sigma_square > 25.)
0218             continue;
0219           if (d0 * d0 / d0_sigma_square < 25.)
0220             multWeighted += constWeight;
0221         } else
0222           multWeighted += constWeight;
0223       } else {  //No track --> neutral particle
0224         if ((constWeight * part->pt()) < 1.0)
0225           continue;  //Only use neutrals with pt > 1 GeV
0226         multWeighted += constWeight;
0227       }
0228     }
0229
0230     float deta = daughter->eta() - jet->eta();
0231     float dphi = reco::deltaPhi(daughter->phi(), jet->phi());
0232     float partPt = constWeight * daughter->pt();
0233     float weight = partPt * partPt;
0234
0235     sum_weight += weight;
0236     sum_pt += partPt;
0237     sum_deta += deta * weight;
0238     sum_dphi += dphi * weight;
0239     sum_deta2 += deta * deta * weight;
0240     sum_detadphi += deta * dphi * weight;
0241     sum_dphi2 += dphi * dphi * weight;
0242   }
0243
0244   //Calculate axis2 and ptD
0245   float a = 0., b = 0., c = 0.;
0246   float ave_deta = 0., ave_dphi = 0., ave_deta2 = 0., ave_dphi2 = 0.;
0247   if (sum_weight > 0) {
0248     ave_deta = sum_deta / sum_weight;
0249     ave_dphi = sum_dphi / sum_weight;
0250     ave_deta2 = sum_deta2 / sum_weight;
0251     ave_dphi2 = sum_dphi2 / sum_weight;
0252     a = ave_deta2 - ave_deta * ave_deta;
0253     b = ave_dphi2 - ave_dphi * ave_dphi;
0254     c = -(sum_detadphi / sum_weight - ave_deta * ave_dphi);
0255   }
0256
0257   int mult = std::round(multWeighted);
0258   float delta = sqrt(fabs((a - b) * (a - b) + 4 * c * c));
0259   float axis2 = (a + b - delta > 0 ? sqrt(0.5 * (a + b - delta)) : 0);
0260   float ptD = (sum_weight > 0 ? sqrt(sum_weight) / sum_pt : 0);
0261   return std::make_tuple(mult, ptD, axis2);
0262 }
0263
0264 /// Descriptions method
0265 void QGTagger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0266   edm::ParameterSetDescription desc;
0267   desc.add<edm::InputTag>("srcJets", edm::InputTag("ak4PFJetsCHS"));
0268   desc.add<edm::InputTag>("srcRho", edm::InputTag("fixedGridRhoFastjetAll"));
0269   desc.add<bool>("computeLikelihood", true);
0270   desc.add<std::string>("jetsLabel", "QGL_AK4PFchs");
0271   desc.add<std::string>("systematicsLabel", "");
0272   desc.add<bool>("useQualityCuts", false);
0273   desc.add<edm::InputTag>("jec", edm::InputTag())->setComment("Jet correction service: only applied when non-empty");
0274   desc.add<edm::InputTag>("srcVertexCollection", edm::InputTag("offlinePrimaryVerticesWithBS"))
0275       ->setComment("Ignored for miniAOD, possible to keep empty");
0276   desc.add<edm::InputTag>("srcConstituentWeights", edm::InputTag())->setComment("Constituent weights ValueMap");
0277   descriptions.add("QGTagger", desc);
0278 }
0279
0280 //define this as a plug-in
0281 DEFINE_FWK_MODULE(QGTagger);