File indexing completed on 2024-04-06 12:24:35
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include <memory>
0018 #include <string>
0019 #include <utility>
0020 #include <cmath>
0021 #include <map>
0022
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/global/EDProducer.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0028 #include "DataFormats/Common/interface/RefToBase.h"
0029 #include "DataFormats/Math/interface/Vector3D.h"
0030 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0031 #include "DataFormats/JetReco/interface/Jet.h"
0032 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0033 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0034 #include "DataFormats/VertexReco/interface/Vertex.h"
0035 #include "DataFormats/Common/interface/ValueMap.h"
0036 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0037 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0038 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0039 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0040 #include "DataFormats/TrackReco/interface/Track.h"
0041 #include "DataFormats/MuonReco/interface/Muon.h"
0042 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0043 #include "DataFormats/BTauReco/interface/SoftLeptonTagInfo.h"
0044 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0045
0046 #include "FWCore/Utilities/interface/InputTag.h"
0047 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0048
0049 #include "FWCore/Framework/interface/ModuleFactory.h"
0050 #include "FWCore/Framework/interface/MakerMacros.h"
0051
0052
0053 #include "DataFormats/Math/interface/LorentzVector.h"
0054 #include "Math/GenVector/PxPyPzM4D.h"
0055 #include "Math/GenVector/VectorUtil.h"
0056 #include "Math/GenVector/Boost.h"
0057
0058 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0059
0060 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0061 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0062 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0063 #include "TrackingTools/IPTools/interface/IPTools.h"
0064
0065 class SoftLepton : public edm::global::EDProducer<> {
0066 public:
0067 explicit SoftLepton(const edm::ParameterSet &iConfig);
0068
0069 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0070
0071 struct TrackCompare {
0072 inline bool operator()(const edm::RefToBase<reco::Track> &t1, const edm::RefToBase<reco::Track> &t2) const {
0073 return t1.key() < t2.key();
0074 }
0075 };
0076
0077 using LeptonIds = std::map<unsigned int, float>;
0078 using Leptons = std::map<edm::RefToBase<reco::Track>, LeptonIds, TrackCompare>;
0079
0080
0081 reco::SoftLeptonTagInfo tag(const edm::RefToBase<reco::Jet> &jet,
0082 const reco::TrackRefVector &tracks,
0083 const Leptons &leptons,
0084 const reco::Vertex &primaryVertex,
0085 const TransientTrackBuilder &builder) const;
0086
0087 protected:
0088
0089
0090 GlobalVector refineJetAxis(const edm::RefToBase<reco::Jet> &jet,
0091 const reco::TrackRefVector &tracks,
0092 const edm::RefToBase<reco::Track> &exclude = edm::RefToBase<reco::Track>()) const;
0093
0094 static double relativeEta(const math::XYZVector &vector, const math::XYZVector &axis);
0095
0096 static double boostedPPar(const math::XYZVector &vector, const math::XYZVector &axis);
0097
0098 private:
0099 void produce(edm::StreamID, edm::Event &event, const edm::EventSetup &setup) const final;
0100
0101
0102 const edm::InputTag m_jets;
0103 const edm::EDGetTokenT<reco::JetTracksAssociationCollection> token_jtas;
0104 const edm::EDGetTokenT<edm::View<reco::Jet> > token_jets;
0105 const edm::EDGetTokenT<reco::VertexCollection> token_primaryVertex;
0106 const edm::InputTag m_leptons;
0107 const edm::EDGetTokenT<edm::View<reco::GsfElectron> > token_gsfElectrons;
0108 const edm::EDGetTokenT<edm::View<reco::Electron> > token_electrons;
0109 const edm::EDGetTokenT<reco::PFCandidateCollection> token_pfElectrons;
0110 const edm::EDGetTokenT<edm::View<reco::Muon> > token_muons;
0111 const edm::EDGetTokenT<edm::View<reco::Track> > token_tracks;
0112 const edm::InputTag m_leptonCands;
0113 const edm::EDGetTokenT<edm::ValueMap<float> > token_leptonCands;
0114 const edm::InputTag m_leptonId;
0115 const edm::EDGetTokenT<edm::ValueMap<float> > token_leptonId;
0116
0117
0118 const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> token_builder;
0119
0120 const edm::EDPutTokenT<reco::SoftLeptonTagInfoCollection> token_put;
0121
0122 const unsigned int m_refineJetAxis;
0123 const double m_deltaRCut;
0124 const double m_chi2Cut;
0125
0126
0127 const muon::SelectionType m_muonSelection;
0128
0129
0130 static const reco::Vertex s_nominalBeamSpot;
0131 };
0132
0133 enum AxisType {
0134 AXIS_CALORIMETRIC = 0,
0135 AXIS_CHARGED_AVERAGE = 1,
0136 AXIS_CHARGED_AVERAGE_NOLEPTON = 2,
0137 AXIS_CHARGED_SUM = 3,
0138 AXIS_CHARGED_SUM_NOLEPTON = 4,
0139 AXIS_CALORIMETRIC_NOLEPTON = 5
0140 };
0141
0142 using namespace std;
0143 using namespace edm;
0144 using namespace reco;
0145 using namespace ROOT::Math::VectorUtil;
0146
0147 typedef edm::View<reco::GsfElectron> GsfElectronView;
0148 typedef edm::View<reco::Electron> ElectronView;
0149 typedef edm::View<reco::Muon> MuonView;
0150
0151
0152 const reco::Vertex SoftLepton::s_nominalBeamSpot(
0153 reco::Vertex::Point(0, 0, 0),
0154 reco::Vertex::Error(ROOT::Math::SVector<double, 6>(0.0015 * 0.0015,
0155 0.0,
0156 0.0015 * 0.0015,
0157 0.0,
0158 0.0,
0159 15. * 15.)),
0160 1,
0161 1,
0162 0);
0163
0164
0165 SoftLepton::SoftLepton(const edm::ParameterSet &iConfig)
0166 : m_jets(iConfig.getParameter<edm::InputTag>("jets")),
0167 token_jtas(mayConsume<reco::JetTracksAssociationCollection>(m_jets)),
0168 token_jets(mayConsume<edm::View<reco::Jet> >(m_jets)),
0169 token_primaryVertex(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("primaryVertex"))),
0170 m_leptons(iConfig.getParameter<edm::InputTag>("leptons")),
0171 token_gsfElectrons(mayConsume<GsfElectronView>(m_leptons)),
0172 token_electrons(mayConsume<ElectronView>(m_leptons)),
0173 token_pfElectrons(mayConsume<reco::PFCandidateCollection>(m_leptons)),
0174 token_muons(mayConsume<MuonView>(m_leptons)),
0175 token_tracks(mayConsume<edm::View<reco::Track> >(m_leptons)),
0176 m_leptonCands(iConfig.getParameter<edm::InputTag>("leptonCands")),
0177 token_leptonCands(mayConsume<edm::ValueMap<float> >(m_leptonCands)),
0178 m_leptonId(iConfig.getParameter<edm::InputTag>("leptonId")),
0179 token_leptonId(mayConsume<edm::ValueMap<float> >(m_leptonId)),
0180 token_builder(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0181 token_put(produces()),
0182 m_refineJetAxis(iConfig.getParameter<unsigned int>("refineJetAxis")),
0183 m_deltaRCut(iConfig.getParameter<double>("leptonDeltaRCut")),
0184 m_chi2Cut(iConfig.getParameter<double>("leptonChi2Cut")),
0185 m_muonSelection((muon::SelectionType)iConfig.getParameter<unsigned int>("muonSelection")) {}
0186
0187
0188 void SoftLepton::produce(edm::StreamID, edm::Event &event, const edm::EventSetup &setup) const {
0189
0190 auto const &transientTrackBuilder = setup.getData(token_builder);
0191
0192
0193
0194
0195 ProductID jets_id;
0196 std::vector<edm::RefToBase<reco::Jet> > jets;
0197 std::vector<reco::TrackRefVector> tracks;
0198 do {
0199 {
0200
0201 edm::Handle<reco::JetTracksAssociationCollection> h_jtas = event.getHandle(token_jtas);
0202 if (h_jtas.isValid()) {
0203 unsigned int size = h_jtas->size();
0204 jets.resize(size);
0205 tracks.resize(size);
0206 for (unsigned int i = 0; i < size; ++i) {
0207 jets[i] = (*h_jtas)[i].first;
0208 tracks[i] = (*h_jtas)[i].second;
0209 }
0210 break;
0211 }
0212 }
0213 {
0214
0215 edm::Handle<edm::View<reco::Jet> > h_jets = event.getHandle(token_jets);
0216 if (h_jets.isValid()) {
0217 unsigned int size = h_jets->size();
0218 jets.resize(size);
0219 tracks.resize(size);
0220 for (unsigned int i = 0; i < h_jets->size(); i++)
0221 jets[i] = h_jets->refAt(i);
0222 break;
0223 }
0224 }
0225 {
0226 throw edm::Exception(edm::errors::NotFound)
0227 << "Object " << m_jets
0228 << " of type among (\"reco::JetTracksAssociationCollection\", \"edm::View<reco::Jet>\") not found";
0229 }
0230 } while (false);
0231
0232
0233 reco::Vertex vertex;
0234 Handle<reco::VertexCollection> h_primaryVertex = event.getHandle(token_primaryVertex);
0235 if (h_primaryVertex.isValid() and not h_primaryVertex->empty())
0236 vertex = h_primaryVertex->front();
0237 else
0238
0239 vertex = s_nominalBeamSpot;
0240
0241
0242 Leptons leptons;
0243
0244 Handle<edm::ValueMap<float> > h_leptonCands;
0245 bool haveLeptonCands = !(m_leptonCands == edm::InputTag());
0246 if (haveLeptonCands)
0247 h_leptonCands = event.getHandle(token_leptonCands);
0248
0249
0250
0251 unsigned int leptonId = SoftLeptonProperties::Quality::leptonId;
0252 do {
0253 {
0254
0255 Handle<GsfElectronView> h_electrons = event.getHandle(token_gsfElectrons);
0256
0257 if (h_electrons.isValid()) {
0258 leptonId = SoftLeptonProperties::Quality::egammaElectronId;
0259 for (GsfElectronView::const_iterator electron = h_electrons->begin(); electron != h_electrons->end();
0260 ++electron) {
0261 LeptonIds &id = leptons[reco::TrackBaseRef(electron->gsfTrack())];
0262 id[SoftLeptonProperties::Quality::pfElectronId] = electron->mva_e_pi();
0263 if (haveLeptonCands)
0264 id[SoftLeptonProperties::Quality::btagElectronCands] =
0265 (*h_leptonCands)[h_electrons->refAt(electron - h_electrons->begin())];
0266 }
0267 break;
0268 }
0269 }
0270 {
0271
0272
0273 Handle<ElectronView> h_electrons = event.getHandle(token_electrons);
0274 if (h_electrons.isValid()) {
0275 leptonId = SoftLeptonProperties::Quality::egammaElectronId;
0276 for (ElectronView::const_iterator electron = h_electrons->begin(); electron != h_electrons->end(); ++electron) {
0277 LeptonIds &id = leptons[reco::TrackBaseRef(electron->track())];
0278 if (haveLeptonCands)
0279 id[SoftLeptonProperties::Quality::btagElectronCands] =
0280 (*h_leptonCands)[h_electrons->refAt(electron - h_electrons->begin())];
0281 }
0282 break;
0283 }
0284 }
0285 {
0286
0287
0288 Handle<reco::PFCandidateCollection> h_electrons = event.getHandle(token_pfElectrons);
0289 if (h_electrons.isValid()) {
0290 leptonId = SoftLeptonProperties::Quality::egammaElectronId;
0291 for (reco::PFCandidateCollection::const_iterator electron = h_electrons->begin();
0292 electron != h_electrons->end();
0293 ++electron) {
0294 LeptonIds *id;
0295 if (electron->gsfTrackRef().isNonnull())
0296 id = &leptons[reco::TrackBaseRef(electron->gsfTrackRef())];
0297 else if (electron->trackRef().isNonnull())
0298 id = &leptons[reco::TrackBaseRef(electron->trackRef())];
0299 else
0300 continue;
0301 (*id)[SoftLeptonProperties::Quality::pfElectronId] = electron->mva_e_pi();
0302 if (haveLeptonCands)
0303 (*id)[SoftLeptonProperties::Quality::btagElectronCands] =
0304 (*h_leptonCands)[reco::PFCandidateRef(h_electrons, electron - h_electrons->begin())];
0305 }
0306 break;
0307 }
0308 }
0309 {
0310
0311 Handle<MuonView> h_muons = event.getHandle(token_muons);
0312 if (h_muons.isValid()) {
0313 for (MuonView::const_iterator muon = h_muons->begin(); muon != h_muons->end(); ++muon) {
0314
0315 if (muon::isGoodMuon(*muon, m_muonSelection)) {
0316 LeptonIds *id;
0317 if (muon->globalTrack().isNonnull())
0318 id = &leptons[reco::TrackBaseRef(muon->globalTrack())];
0319 else if (muon->innerTrack().isNonnull())
0320 id = &leptons[reco::TrackBaseRef(muon->innerTrack())];
0321 else if (muon->outerTrack().isNonnull())
0322
0323 id = &leptons[reco::TrackBaseRef(muon->outerTrack())];
0324 else
0325 continue;
0326 if (haveLeptonCands)
0327 (*id)[SoftLeptonProperties::Quality::btagMuonCands] =
0328 (*h_leptonCands)[h_muons->refAt(muon - h_muons->begin())];
0329 }
0330 }
0331 break;
0332 }
0333 }
0334 {
0335
0336 Handle<edm::View<reco::Track> > h_tracks = event.getHandle(token_tracks);
0337 if (h_tracks.isValid()) {
0338 for (unsigned int i = 0; i < h_tracks->size(); i++) {
0339 LeptonIds &id = leptons[h_tracks->refAt(i)];
0340 if (haveLeptonCands)
0341 id[SoftLeptonProperties::Quality::btagLeptonCands] = (*h_leptonCands)[h_tracks->refAt(i)];
0342 }
0343 break;
0344 }
0345 }
0346 {
0347 throw edm::Exception(edm::errors::NotFound) << "Object " << m_leptons
0348 << " of type among (\"edm::View<reco::GsfElectron>\", "
0349 "\"edm::View<reco::Muon>\", \"edm::View<reco::Track>\") !found";
0350 }
0351 } while (false);
0352
0353 if (!(m_leptonId == edm::InputTag())) {
0354 edm::ValueMap<float> const &h_leptonId = event.get(token_leptonId);
0355
0356 for (Leptons::iterator lepton = leptons.begin(); lepton != leptons.end(); ++lepton)
0357 lepton->second[leptonId] = h_leptonId[lepton->first];
0358 }
0359
0360
0361 reco::SoftLeptonTagInfoCollection outputCollection;
0362 for (unsigned int i = 0; i < jets.size(); ++i) {
0363 reco::SoftLeptonTagInfo result = tag(jets[i], tracks[i], leptons, vertex, transientTrackBuilder);
0364 outputCollection.push_back(result);
0365 }
0366 event.emplace(token_put, std::move(outputCollection));
0367 }
0368
0369
0370 reco::SoftLeptonTagInfo SoftLepton::tag(const edm::RefToBase<reco::Jet> &jet,
0371 const reco::TrackRefVector &tracks,
0372 const Leptons &leptons,
0373 const reco::Vertex &primaryVertex,
0374 const TransientTrackBuilder &transientTrackBuilder) const {
0375 reco::SoftLeptonTagInfo info;
0376 info.setJetRef(jet);
0377
0378 for (Leptons::const_iterator lepton = leptons.begin(); lepton != leptons.end(); ++lepton) {
0379 const math::XYZVector &lepton_momentum = lepton->first->momentum();
0380 if (m_chi2Cut > 0.0 && lepton->first->normalizedChi2() > m_chi2Cut)
0381 continue;
0382
0383 const GlobalVector jetAxis = refineJetAxis(jet, tracks, lepton->first);
0384 const math::XYZVector axis(jetAxis.x(), jetAxis.y(), jetAxis.z());
0385 float deltaR = Geom::deltaR(lepton_momentum, axis);
0386 if (deltaR > m_deltaRCut)
0387 continue;
0388
0389 reco::SoftLeptonProperties properties;
0390
0391 reco::TransientTrack transientTrack = transientTrackBuilder.build(*lepton->first);
0392 Measurement1D ip2d = IPTools::signedTransverseImpactParameter(transientTrack, jetAxis, primaryVertex).second;
0393 Measurement1D ip3d = IPTools::signedImpactParameter3D(transientTrack, jetAxis, primaryVertex).second;
0394 properties.sip2dsig = ip2d.significance();
0395 properties.sip3dsig = ip3d.significance();
0396 properties.sip2d = ip2d.value();
0397 properties.sip3d = ip3d.value();
0398 properties.deltaR = deltaR;
0399 properties.ptRel = Perp(lepton_momentum, axis);
0400 properties.p0Par = boostedPPar(lepton_momentum, axis);
0401 properties.etaRel = relativeEta(lepton_momentum, axis);
0402 properties.ratio = lepton_momentum.R() / axis.R();
0403 properties.ratioRel = lepton_momentum.Dot(axis) / axis.Mag2();
0404
0405 for (LeptonIds::const_iterator iter = lepton->second.begin(); iter != lepton->second.end(); ++iter)
0406 properties.setQuality(static_cast<SoftLeptonProperties::Quality::Generic>(iter->first), iter->second);
0407
0408 info.insert(lepton->first, properties);
0409 }
0410
0411 return info;
0412 }
0413
0414
0415 GlobalVector SoftLepton::refineJetAxis(const edm::RefToBase<reco::Jet> &jet,
0416 const reco::TrackRefVector &tracks,
0417 const reco::TrackBaseRef &exclude
0418 ) const {
0419 math::XYZVector axis = jet->momentum();
0420
0421 if (m_refineJetAxis == AXIS_CHARGED_AVERAGE or m_refineJetAxis == AXIS_CHARGED_AVERAGE_NOLEPTON) {
0422 double sum_pT = 0.;
0423 double sum_eta_by_pT = 0.;
0424 double sum_phi_by_pT = 0.;
0425
0426 double perp;
0427 double phi_rel;
0428 double eta_rel;
0429
0430
0431 for (reco::TrackRefVector::const_iterator track_it = tracks.begin(); track_it != tracks.end(); ++track_it) {
0432 const reco::Track &track = **track_it;
0433
0434 perp = track.pt();
0435 eta_rel = (double)track.eta() - axis.eta();
0436 phi_rel = (double)track.phi() - axis.phi();
0437 while (phi_rel < -M_PI)
0438 phi_rel += 2 * M_PI;
0439 while (phi_rel > M_PI)
0440 phi_rel -= 2 * M_PI;
0441
0442 sum_pT += perp;
0443 sum_phi_by_pT += perp * phi_rel;
0444 sum_eta_by_pT += perp * eta_rel;
0445 }
0446
0447
0448 if (m_refineJetAxis == AXIS_CHARGED_AVERAGE_NOLEPTON and exclude.isNonnull()) {
0449 const reco::Track &track = *exclude;
0450
0451 perp = track.pt();
0452 eta_rel = (double)track.eta() - axis.eta();
0453 phi_rel = (double)track.phi() - axis.phi();
0454 while (phi_rel < -M_PI)
0455 phi_rel += 2 * M_PI;
0456 while (phi_rel > M_PI)
0457 phi_rel -= 2 * M_PI;
0458
0459 sum_pT -= perp;
0460 sum_phi_by_pT -= perp * phi_rel;
0461 sum_eta_by_pT -= perp * eta_rel;
0462 }
0463
0464 if (sum_pT > 1.)
0465 axis =
0466 math::RhoEtaPhiVector(axis.rho(), axis.eta() + sum_eta_by_pT / sum_pT, axis.phi() + sum_phi_by_pT / sum_pT);
0467
0468 } else if (m_refineJetAxis == AXIS_CHARGED_SUM or m_refineJetAxis == AXIS_CHARGED_SUM_NOLEPTON) {
0469 math::XYZVector sum;
0470
0471
0472 for (reco::TrackRefVector::const_iterator track_it = tracks.begin(); track_it != tracks.end(); ++track_it) {
0473 const reco::Track &track = **track_it;
0474 sum += track.momentum();
0475 }
0476
0477
0478 if (m_refineJetAxis == AXIS_CHARGED_SUM_NOLEPTON and exclude.isNonnull()) {
0479 const reco::Track &track = *exclude;
0480 sum -= track.momentum();
0481 }
0482
0483 if (sum.R() > 1.)
0484 axis = sum;
0485 } else if (m_refineJetAxis == AXIS_CALORIMETRIC_NOLEPTON) {
0486 axis -= exclude->momentum();
0487 }
0488
0489 return GlobalVector(axis.x(), axis.y(), axis.z());
0490 }
0491
0492 double SoftLepton::relativeEta(const math::XYZVector &vector, const math::XYZVector &axis) {
0493 double mag = vector.r() * axis.r();
0494 double dot = vector.Dot(axis);
0495 return -log((mag - dot) / (mag + dot)) / 2;
0496 }
0497
0498
0499 double SoftLepton::boostedPPar(const math::XYZVector &vector, const math::XYZVector &axis) {
0500 static const double lepton_mass = 0.00;
0501 static const double jet_mass = 5.279;
0502 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > lepton(
0503 vector.Dot(axis) / axis.r(), Perp(vector, axis), 0., lepton_mass);
0504 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double> > jet(axis.r(), 0., 0., jet_mass);
0505 ROOT::Math::BoostX boost(-jet.Beta());
0506 return boost(lepton).x();
0507 }
0508
0509
0510 void SoftLepton::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0511 edm::ParameterSetDescription desc;
0512 desc.add<unsigned int>("muonSelection", 1);
0513 desc.add<edm::InputTag>("leptons", edm::InputTag("muons"));
0514 desc.add<edm::InputTag>("primaryVertex", edm::InputTag("offlinePrimaryVertices"));
0515 desc.add<edm::InputTag>("leptonCands", edm::InputTag());
0516 desc.add<edm::InputTag>("leptonId", edm::InputTag());
0517 desc.add<unsigned int>("refineJetAxis", 0);
0518 desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
0519 desc.add<double>("leptonDeltaRCut", 0.4);
0520 desc.add<double>("leptonChi2Cut", 9999.0);
0521 descriptions.addDefault(desc);
0522 }
0523
0524 DEFINE_FWK_MODULE(SoftLepton);