Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-23 03:25:23

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       paramsToken(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("jetsLabel")))),
0031       useQC(iConfig.getParameter<bool>("useQualityCuts")),
0032       useJetCorr(!iConfig.getParameter<edm::InputTag>("jec").label().empty()),
0033       produceSyst(!iConfig.getParameter<std::string>("systematicsLabel").empty()) {
0034   produces<edm::ValueMap<float>>("qgLikelihood");
0035   produces<edm::ValueMap<float>>("axis2");
0036   produces<edm::ValueMap<int>>("mult");
0037   produces<edm::ValueMap<float>>("ptD");
0038   if (produceSyst) {
0039     systToken = esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("systematicsLabel")));
0040     produces<edm::ValueMap<float>>("qgLikelihoodSmearedQuark");
0041     produces<edm::ValueMap<float>>("qgLikelihoodSmearedGluon");
0042     produces<edm::ValueMap<float>>("qgLikelihoodSmearedAll");
0043   }
0044 }
0045 
0046 /// Produce qgLikelihood using {mult, ptD, -log(axis2)}
0047 void QGTagger::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0048   std::vector<float> qgProduct;
0049   std::vector<float> axis2Product;
0050   std::vector<int> multProduct;
0051   std::vector<float> ptDProduct;
0052   std::vector<float> smearedQuarkProduct;
0053   std::vector<float> smearedGluonProduct;
0054   std::vector<float> smearedAllProduct;
0055 
0056   edm::Handle<edm::View<reco::Jet>> jets;
0057   iEvent.getByToken(jetsToken, jets);
0058   edm::Handle<reco::JetCorrector> jetCorr;
0059   if (useJetCorr)
0060     iEvent.getByToken(jetCorrectorToken, jetCorr);
0061   edm::Handle<reco::VertexCollection> vertexCollection;
0062   iEvent.getByToken(vertexToken, vertexCollection);
0063   edm::Handle<double> rho;
0064   iEvent.getByToken(rhoToken, rho);
0065 
0066   const auto& QGLParamsColl = iSetup.getData(paramsToken);
0067 
0068   const QGLikelihoodSystematicsObject* QGLSystColl = nullptr;
0069   if (produceSyst) {
0070     QGLSystColl = &iSetup.getData(systToken);
0071   }
0072 
0073   bool weStillNeedToCheckJetCandidates = true;
0074   bool weAreUsingPackedCandidates = false;
0075   for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
0076     if (weStillNeedToCheckJetCandidates) {
0077       weAreUsingPackedCandidates = isPackedCandidate(&*jet);
0078       weStillNeedToCheckJetCandidates = false;
0079     }
0080     float pt = (useJetCorr ? jet->pt() * jetCorr->correction(*jet) : jet->pt());
0081 
0082     float ptD, axis2;
0083     int mult;
0084     std::tie(mult, ptD, axis2) = calcVariables(&*jet, vertexCollection, weAreUsingPackedCandidates);
0085 
0086     float qgValue;
0087     if (mult > 2)
0088       qgValue =
0089           qgLikelihood.computeQGLikelihood(QGLParamsColl, pt, jet->eta(), *rho, {(float)mult, ptD, -std::log(axis2)});
0090     else
0091       qgValue = -1;
0092 
0093     qgProduct.push_back(qgValue);
0094     if (produceSyst) {
0095       smearedQuarkProduct.push_back(qgLikelihood.systematicSmearing(*QGLSystColl, pt, jet->eta(), *rho, qgValue, 0));
0096       smearedGluonProduct.push_back(qgLikelihood.systematicSmearing(*QGLSystColl, pt, jet->eta(), *rho, qgValue, 1));
0097       smearedAllProduct.push_back(qgLikelihood.systematicSmearing(*QGLSystColl, pt, jet->eta(), *rho, qgValue, 2));
0098     }
0099     axis2Product.push_back(axis2);
0100     multProduct.push_back(mult);
0101     ptDProduct.push_back(ptD);
0102   }
0103 
0104   putInEvent("qgLikelihood", jets, qgProduct, iEvent);
0105   putInEvent("axis2", jets, axis2Product, iEvent);
0106   putInEvent("mult", jets, multProduct, iEvent);
0107   putInEvent("ptD", jets, ptDProduct, iEvent);
0108   if (produceSyst) {
0109     putInEvent("qgLikelihoodSmearedQuark", jets, smearedQuarkProduct, iEvent);
0110     putInEvent("qgLikelihoodSmearedGluon", jets, smearedGluonProduct, iEvent);
0111     putInEvent("qgLikelihoodSmearedAll", jets, smearedAllProduct, iEvent);
0112   }
0113 }
0114 
0115 /// Function to put product into event
0116 template <typename T>
0117 void QGTagger::putInEvent(const std::string& name,
0118                           const edm::Handle<edm::View<reco::Jet>>& jets,
0119                           const std::vector<T>& product,
0120                           edm::Event& iEvent) const {
0121   auto out = std::make_unique<edm::ValueMap<T>>();
0122   typename edm::ValueMap<T>::Filler filler(*out);
0123   filler.insert(jets, product.begin(), product.end());
0124   filler.fill();
0125   iEvent.put(std::move(out), name);
0126 }
0127 
0128 /// Function to tell us if we are using packedCandidates, only test for first candidate
0129 bool QGTagger::isPackedCandidate(const reco::Jet* jet) const {
0130   for (auto candidate : jet->getJetConstituentsQuick()) {
0131     if (typeid(pat::PackedCandidate) == typeid(*candidate))
0132       return true;
0133     else if (typeid(reco::PFCandidate) == typeid(*candidate))
0134       return false;
0135     else
0136       throw cms::Exception("WrongJetCollection", "Jet constituents are not particle flow candidates");
0137   }
0138   return false;
0139 }
0140 
0141 /// Calculation of axis2, mult and ptD
0142 std::tuple<int, float, float> QGTagger::calcVariables(const reco::Jet* jet,
0143                                                       edm::Handle<reco::VertexCollection>& vC,
0144                                                       bool weAreUsingPackedCandidates) const {
0145   float sum_weight = 0., sum_deta = 0., sum_dphi = 0., sum_deta2 = 0., sum_dphi2 = 0., sum_detadphi = 0., sum_pt = 0.;
0146   int mult = 0;
0147 
0148   //Loop over the jet constituents
0149   for (auto daughter : jet->getJetConstituentsQuick()) {
0150     if (weAreUsingPackedCandidates) {  //packed candidate situation
0151       auto part = static_cast<const pat::PackedCandidate*>(daughter);
0152 
0153       if (part->charge()) {
0154         if (!(part->fromPV() > 1 && part->trackHighPurity()))
0155           continue;
0156         if (useQC) {
0157           if ((part->dz() * part->dz()) / (part->dzError() * part->dzError()) > 25.)
0158             continue;
0159           if ((part->dxy() * part->dxy()) / (part->dxyError() * part->dxyError()) < 25.)
0160             ++mult;
0161         } else
0162           ++mult;
0163       } else {
0164         if (part->pt() < 1.0)
0165           continue;
0166         ++mult;
0167       }
0168     } else {
0169       auto part = static_cast<const reco::PFCandidate*>(daughter);
0170 
0171       reco::TrackRef itrk = part->trackRef();
0172       if (itrk.isNonnull()) {  //Track exists --> charged particle
0173         auto vtxLead = vC->begin();
0174         auto vtxClose = vC->begin();  //Search for closest vertex to track
0175         for (auto vtx = vC->begin(); vtx != vC->end(); ++vtx) {
0176           if (fabs(itrk->dz(vtx->position())) < fabs(itrk->dz(vtxClose->position())))
0177             vtxClose = vtx;
0178         }
0179         if (!(vtxClose == vtxLead && itrk->quality(reco::TrackBase::qualityByName("highPurity"))))
0180           continue;
0181 
0182         if (useQC) {  //If useQC, require dz and d0 cuts
0183           float dz = itrk->dz(vtxClose->position());
0184           float d0 = itrk->dxy(vtxClose->position());
0185           float dz_sigma_square = pow(itrk->dzError(), 2) + pow(vtxClose->zError(), 2);
0186           float d0_sigma_square = pow(itrk->d0Error(), 2) + pow(vtxClose->xError(), 2) + pow(vtxClose->yError(), 2);
0187           if (dz * dz / dz_sigma_square > 25.)
0188             continue;
0189           if (d0 * d0 / d0_sigma_square < 25.)
0190             ++mult;
0191         } else
0192           ++mult;
0193       } else {  //No track --> neutral particle
0194         if (part->pt() < 1.0)
0195           continue;  //Only use neutrals with pt > 1 GeV
0196         ++mult;
0197       }
0198     }
0199 
0200     float deta = daughter->eta() - jet->eta();
0201     float dphi = reco::deltaPhi(daughter->phi(), jet->phi());
0202     float partPt = daughter->pt();
0203     float weight = partPt * partPt;
0204 
0205     sum_weight += weight;
0206     sum_pt += partPt;
0207     sum_deta += deta * weight;
0208     sum_dphi += dphi * weight;
0209     sum_deta2 += deta * deta * weight;
0210     sum_detadphi += deta * dphi * weight;
0211     sum_dphi2 += dphi * dphi * weight;
0212   }
0213 
0214   //Calculate axis2 and ptD
0215   float a = 0., b = 0., c = 0.;
0216   float ave_deta = 0., ave_dphi = 0., ave_deta2 = 0., ave_dphi2 = 0.;
0217   if (sum_weight > 0) {
0218     ave_deta = sum_deta / sum_weight;
0219     ave_dphi = sum_dphi / sum_weight;
0220     ave_deta2 = sum_deta2 / sum_weight;
0221     ave_dphi2 = sum_dphi2 / sum_weight;
0222     a = ave_deta2 - ave_deta * ave_deta;
0223     b = ave_dphi2 - ave_dphi * ave_dphi;
0224     c = -(sum_detadphi / sum_weight - ave_deta * ave_dphi);
0225   }
0226   float delta = sqrt(fabs((a - b) * (a - b) + 4 * c * c));
0227   float axis2 = (a + b - delta > 0 ? sqrt(0.5 * (a + b - delta)) : 0);
0228   float ptD = (sum_weight > 0 ? sqrt(sum_weight) / sum_pt : 0);
0229   return std::make_tuple(mult, ptD, axis2);
0230 }
0231 
0232 /// Descriptions method
0233 void QGTagger::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0234   edm::ParameterSetDescription desc;
0235   desc.add<edm::InputTag>("srcJets");
0236   desc.add<edm::InputTag>("srcRho");
0237   desc.add<std::string>("jetsLabel");
0238   desc.add<std::string>("systematicsLabel", "");
0239   desc.add<bool>("useQualityCuts");
0240   desc.add<edm::InputTag>("jec", edm::InputTag())->setComment("Jet correction service: only applied when non-empty");
0241   desc.add<edm::InputTag>("srcVertexCollection")->setComment("Ignored for miniAOD, possible to keep empty");
0242   descriptions.add("QGTagger", desc);
0243 }
0244 
0245 //define this as a plug-in
0246 DEFINE_FWK_MODULE(QGTagger);