From 8901e721728c21b112d5c01e207adb7d4590f88b Mon Sep 17 00:00:00 2001 From: Lovisa Date: Mon, 8 Dec 2025 14:46:53 +0100 Subject: [PATCH 1/7] Including EXOnanoAODv1 custom constent --- .../PyReleaseValidation/python/relval_nano.py | 12 + .../NanoAOD/interface/DisTauTagScaling.h | 191 ++++ PhysicsTools/NanoAOD/plugins/BuildFile.xml | 10 + .../NanoAOD/plugins/DSAMuonTableProducer.cc | 318 +++++++ PhysicsTools/NanoAOD/plugins/DisTauTag.cc | 338 +++++++ .../NanoAOD/plugins/DispJetTableProducer.cc | 491 ++++++++++ .../plugins/ElectronExtendedTableProducer.cc | 295 ++++++ .../plugins/ElectronVertexTableProducer.cc | 468 ++++++++++ .../NanoAOD/plugins/JetImpactParameters.cc | 140 +++ .../plugins/MuonExtendedTableProducer.cc | 338 +++++++ .../plugins/MuonVertexTableProducer.cc | 874 ++++++++++++++++++ .../plugins/cscMDSshowerTableProducer.cc | 390 ++++++++ .../plugins/dtMDSshowerTableProducer.cc | 339 +++++++ PhysicsTools/NanoAOD/python/autoNANO.py | 3 + .../NanoAOD/python/custom_displacedtau_cff.py | 157 ++++ PhysicsTools/NanoAOD/python/custom_exo_cff.py | 315 +++++++ PhysicsTools/NanoAOD/test/test-exoNano.sh | 1 + 17 files changed, 4680 insertions(+) create mode 100644 PhysicsTools/NanoAOD/interface/DisTauTagScaling.h create mode 100644 PhysicsTools/NanoAOD/plugins/DSAMuonTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/plugins/DisTauTag.cc create mode 100644 PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/plugins/ElectronVertexTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/plugins/JetImpactParameters.cc create mode 100644 PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/plugins/MuonVertexTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/plugins/dtMDSshowerTableProducer.cc create mode 100644 PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py create mode 100644 PhysicsTools/NanoAOD/python/custom_exo_cff.py create mode 100755 PhysicsTools/NanoAOD/test/test-exoNano.sh diff --git a/Configuration/PyReleaseValidation/python/relval_nano.py b/Configuration/PyReleaseValidation/python/relval_nano.py index 476ab6a6333c8..c25931b6d2a12 100644 --- a/Configuration/PyReleaseValidation/python/relval_nano.py +++ b/Configuration/PyReleaseValidation/python/relval_nano.py @@ -306,6 +306,9 @@ def subnext(self): steps['BTVNANO_mc15.0'] = merge([{'-s': 'NANO:@BTV,DQM:@nanoAODDQM', '-n': '1000'}, steps['NANO_mc15.0']]) +steps['EXONANO_mc15.0'] = merge([{'-s': 'PAT,NANO:@EXO', '-n': '1000'}, + steps['NANO_mc15.0']]) + steps['lepTrackInfoNANO_mc15.0'] = merge([{'-s': 'NANO:@LepTrackInfo,DQM:@nanoAODDQM', '-n': '1000'}, steps['NANO_mc15.0']]) @@ -323,10 +326,14 @@ def subnext(self): # ====== DATA ====== lumis_Run2025C = {392175: [[95, 542]]} +lumis_Run2025C_Muon0 = {392997: [[146, 202]]} steps['JetMET1_Run2025C_MINIAOD_150X'] = {'INPUT': InputInfo( location='STD', ls=lumis_Run2025C, dataSet='/JetMET1/Run2025C-PromptReco-v1/MINIAOD')} +steps['Muon0_Run2025C_AOD_150X'] = {'INPUT': InputInfo( + location='STD',ls=lumis_Run2025C_Muon0, dataSet='/Muon0/Run2025C-PromptReco-v1/AOD')} + steps['ScoutingPFRun3_Run2025C_HLTSCOUT_150X'] = {'INPUT': InputInfo(location='STD', ls=lumis_Run2025C, dataSet='/ScoutingPFRun3/Run2025C-v1/HLTSCOUT')} @@ -350,6 +357,9 @@ def subnext(self): steps['BTVNANO_data15.0'] = merge([{'-s': 'NANO:@BTV,DQM:@nanoAODDQM', '-n': '1000'}, steps['NANO_data15.0']]) +steps['EXONANO_data15.0'] = merge([{'-s': 'PAT,NANO:@EXO', '-n': '1000'}, + steps['NANO_data15.0']]) + steps['lepTrackInfoNANO_data15.0'] = merge([{'-s': 'NANO:@LepTrackInfo,DQM:@nanoAODDQM', '-n': '1000'}, steps['NANO_data15.0']]) @@ -495,6 +505,7 @@ def subnext(self): workflows[_wfn()] = ['ScoutingNANOmc150X', ['TTbar_13p6_Summer24_MINIAOD_150X', 'scoutingNANO_mc15.0']] workflows[_wfn()] = ['ScoutingNANOwithPromptmc150X', ['TTbar_13p6_Summer24_MINIAOD_150X', 'scoutingNANO_withPrompt_mc15.0']] workflows[_wfn()] = ['BPHNANOmc150X', ['TTbar_13p6_Summer24_MINIAOD_150X', 'BPHNANO_mc15.0']] +workflows[_wfn()] = ['EXONANOmc150X', ['TTbar_13p6_Summer24_AOD_140X', 'EXONANO_mc15.0']] # POG/PAG custom NANOs, data _wfn.subnext() @@ -507,6 +518,7 @@ def subnext(self): workflows[_wfn()] = ['ScoutingNANOdata150Xrun3', ['ScoutingPFRun3_Run2025C_HLTSCOUT_150X', 'scoutingNANO_data15.0']] workflows[_wfn()] = ['ScoutingNANOwithPromptdata150Xrun3', ['ScoutingPFMonitor_Run2025C_MINIAOD_150X', 'scoutingNANO_withPrompt_data15.0']] # noqa workflows[_wfn()] = ['BPHNANOdata150Xrun3', ['JetMET1_Run2025C_MINIAOD_150X', 'BPHNANO_data15.0']] +workflows[_wfn()] = ['EXONANOdata150Xrun3', ['Muon0_Run2025C_AOD_150X', 'EXONANO_data15.0']] # DPG custom NANOs, data _wfn.subnext() diff --git a/PhysicsTools/NanoAOD/interface/DisTauTagScaling.h b/PhysicsTools/NanoAOD/interface/DisTauTagScaling.h new file mode 100644 index 0000000000000..8027dc6a679a9 --- /dev/null +++ b/PhysicsTools/NanoAOD/interface/DisTauTagScaling.h @@ -0,0 +1,191 @@ +namespace scaling { + struct PfCand { + inline static const std::vector> mean = {{0, 0}, + {50.0, 50.0}, + {0.0, 0.0}, + {0.0, 0.0}, + {0.09469, 0.09469}, + {0, 0}, + {0, 0}, + {0, 0}, + {5.0, 5.0}, + {5.0, 5.0}, + {6.442, 6.442}, + {0, 0}, + {-0.01635, -0.01635}, + {0.02704, 0.02704}, + {-0.01443, -0.01443}, + {0.05551, 0.05551}, + {11.16, 11.16}, + {13.53, 13.53}, + {0.5, 0.5}, + {0.5, 0.5}, + {1.3, 1.3}, + {0.5, 0.5}, + {0, 0}, + {0, 0}}; + inline static const std::vector> std = {{1, 1}, {50.0, 50.0}, + {3.0, 3.0}, {3.141592653589793, 3.141592653589793}, + {0.0651, 0.0651}, {1, 1}, + {1, 1}, {1, 1}, + {5.0, 5.0}, {5.0, 5.0}, + {8.344, 8.344}, {1, 1}, + {2.4, 2.4}, {0.08913, 0.08913}, + {7.444, 7.444}, {0.5998, 0.5998}, + {60.39, 60.39}, {6.44, 6.44}, + {0.5, 0.5}, {0.5, 0.5}, + {1.3, 1.3}, {0.5, 0.5}, + {1, 1}, {1, 1}}; + inline static const std::vector> lim_min = { + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-1.0, -1.0}, + {-1.0, -1.0}, + {-1.0, -1.0}, + {-5, -5}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-1.0, -1.0}, + {-1.0, -1.0}, + {-5, -5}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-5, -5}, + {-5, -5}, + {-5, -5}, + {-5, -5}, + {-5, -5}, + {-5, -5}, + {-1.0, -1.0}, + {-1.0, -1.0}, + {-1.0, -1.0}, + {-1.0, -1.0}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}}; + inline static const std::vector> lim_max = { + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {1.0, 1.0}, + {1.0, 1.0}, + {1.0, 1.0}, + {5, 5}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {1.0, 1.0}, + {1.0, 1.0}, + {5, 5}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {5, 5}, + {5, 5}, + {5, 5}, + {5, 5}, + {5, 5}, + {5, 5}, + {1.0, 1.0}, + {1.0, 1.0}, + {1.0, 1.0}, + {1.0, 1.0}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}}; + }; + struct PfCandCategorical { + inline static const std::vector> mean = {{0, 0}, {0, 0}, {0, 0}}; + inline static const std::vector> std = {{1, 1}, {1, 1}, {1, 1}}; + inline static const std::vector> lim_min = { + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}, + {-std::numeric_limits::infinity(), -std::numeric_limits::infinity()}}; + inline static const std::vector> lim_max = { + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}, + {std::numeric_limits::infinity(), std::numeric_limits::infinity()}}; + }; +}; // namespace scaling + +namespace setupSizes { + const inline long int output_classes = 2; + const inline size_t n_Global = 5; + const inline size_t n_PfCand = 24; + const inline size_t nSeq_PfCand = 50; + const inline size_t n_PfCandCategorical = 3; + const inline size_t nSeq_PfCandCategorical = 50; + const inline std::vector CellObjectTypes{"PfCand", "PfCandCategorical"}; +}; // namespace setupSizes + +enum class Global_Features { jet_pt = 0, jet_eta = 1, Lxy = 2, Lz = 3, Lrel = 4 }; + +enum class PfCand_Features { + pfCand_valid = 0, + pfCand_pt = 1, + pfCand_eta = 2, + pfCand_phi = 3, + pfCand_mass = 4, + pfCand_charge = 5, + pfCand_puppiWeight = 6, + pfCand_puppiWeightNoLep = 7, + pfCand_lostInnerHits = 8, + pfCand_nPixelHits = 9, + pfCand_nHits = 10, + pfCand_hasTrackDetails = 11, + pfCand_dxy = 12, + pfCand_dxy_error = 13, + pfCand_dz = 14, + pfCand_dz_error = 15, + pfCand_track_chi2 = 16, + pfCand_track_ndof = 17, + pfCand_caloFraction = 18, + pfCand_hcalFraction = 19, + pfCand_rawCaloFraction = 20, + pfCand_rawHcalFraction = 21, + pfCand_deta = 22, + pfCand_dphi = 23 +}; + +enum class PfCandCategorical_Features { pfCand_particleType = 0, pfCand_pvAssociationQuality = 1, pfCand_fromPV = 2 }; + +enum class CellObjectType { PfCand, PfCandCategorical }; + +template +struct FeaturesHelper; + +template <> +struct FeaturesHelper { + static constexpr CellObjectType object_type = CellObjectType::PfCand; + static constexpr size_t size = 24; + static constexpr size_t length = 50; + using scaler_type = scaling::PfCand; +}; + +template <> +struct FeaturesHelper { + static constexpr CellObjectType object_type = CellObjectType::PfCandCategorical; + static constexpr size_t size = 3; + static constexpr size_t length = 50; + using scaler_type = scaling::PfCandCategorical; +}; + +using FeatureTuple = std::tuple; + +enum class PFParticleType { + Undefined = 0, // undefined + h = 1, // charged hadron + e = 2, // electron + mu = 3, // muon + gamma = 4, // photon + h0 = 5, // neutral hadron + h_HF = 6, // HF tower identified as a hadron + egamma_HF = 7, // HF tower identified as an EM particle +}; + +inline PFParticleType TranslatePdgIdToPFParticleType(int pdgId) { + static const std::map type_map = { + {11, PFParticleType::e}, + {13, PFParticleType::mu}, + {22, PFParticleType::gamma}, + {211, PFParticleType::h}, + {130, PFParticleType::h0}, + {1, PFParticleType::h_HF}, + {2, PFParticleType::egamma_HF}, + }; + auto iter = type_map.find(std::abs(pdgId)); + return iter == type_map.end() ? PFParticleType::Undefined : iter->second; +} diff --git a/PhysicsTools/NanoAOD/plugins/BuildFile.xml b/PhysicsTools/NanoAOD/plugins/BuildFile.xml index bbe8818cfa63a..003e5b6107791 100644 --- a/PhysicsTools/NanoAOD/plugins/BuildFile.xml +++ b/PhysicsTools/NanoAOD/plugins/BuildFile.xml @@ -22,8 +22,17 @@ + + + + + + + + + @@ -34,6 +43,7 @@ + diff --git a/PhysicsTools/NanoAOD/plugins/DSAMuonTableProducer.cc b/PhysicsTools/NanoAOD/plugins/DSAMuonTableProducer.cc new file mode 100644 index 0000000000000..974522e02a33d --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/DSAMuonTableProducer.cc @@ -0,0 +1,318 @@ +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonSelectors.h" +#include "DataFormats/NanoAOD/interface/FlatTable.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/IPTools/interface/IPTools.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +class DSAMuonTableProducer : public edm::global::EDProducer<> { +public: + DSAMuonTableProducer(const edm::ParameterSet& iConfig) + : name_(iConfig.getParameter("name")), + dsaMuonTag_(consumes>(iConfig.getParameter("dsaMuons"))), + muonTag_(consumes>(iConfig.getParameter("muons"))), + vtxTag_(consumes(iConfig.getParameter("primaryVertex"))), + bsTag_(consumes(iConfig.getParameter("beamspot"))), + transientTrackBuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) { + produces(); + } + + ~DSAMuonTableProducer() override {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("name")->setComment("name of the muon nanoaod::FlatTable"); + desc.add("dsaMuons")->setComment("input displaced standalone muon collection"); + desc.add("muons")->setComment("input muon collection"); + desc.add("primaryVertex")->setComment("input primary vertex collection"); + desc.add("beamspot")->setComment("input beamspot collection"); + descriptions.add("DSAMuonTable", desc); + } + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + bool passesDisplacedID(const reco::Track& dsaMuon) const; + std::tuple getMatches(const pat::Muon& muon, + const reco::Track& dsaMuon, + float minPositionDiff = 1e-6) const; + std::tuple getSegments(const reco::Track& dsaMuon) const; + + std::string name_; + edm::EDGetTokenT> dsaMuonTag_; + edm::EDGetTokenT> muonTag_; + edm::EDGetTokenT vtxTag_; + edm::EDGetTokenT bsTag_; + edm::ESGetToken transientTrackBuilderToken_; +}; + +void DSAMuonTableProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + float minPositionDiffForMatching = 1e-6; + + const std::vector& dsaMuons = iEvent.get(dsaMuonTag_); + const std::vector& patMuons = iEvent.get(muonTag_); + + const reco::VertexCollection& primaryVertices = iEvent.get(vtxTag_); + const auto& pv = primaryVertices.at(0); + GlobalPoint primaryVertex(pv.x(), pv.y(), pv.z()); + + const reco::BeamSpot& beamSpotInput = iEvent.get(bsTag_); + const auto& bs = beamSpotInput.position(); + GlobalPoint beamSpot(bs.x(), bs.y(), bs.z()); + reco::Vertex beamSpotVertex(beamSpotInput.position(), beamSpotInput.covariance3D()); + + const TransientTrackBuilder& builder = iSetup.getData(transientTrackBuilderToken_); + + unsigned int nDSAMuons = dsaMuons.size(); + unsigned int nMuons = patMuons.size(); + + std::vector idx, trkNumPlanes, trkNumHits, trkNumDTHits, trkNumCSCHits, outerEta, outerPhi; + + std::vector dzPV, dzPVErr, dxyPVTraj, dxyPVTrajErr, dxyPVSigned, dxyPVSignedErr, ip3DPVSigned, ip3DPVSignedErr; + + std::vector displacedId; + std::vector> nMatchesPerMuon; + std::vector muonMatch1, muonMatch1idx, muonMatch2, muonMatch2idx, muonMatch3, muonMatch3idx, muonMatch4, + muonMatch4idx, muonMatch5, muonMatch5idx; + std::vector muonDTMatch1, muonDTMatch1idx, muonDTMatch2, muonDTMatch2idx, muonDTMatch3, muonDTMatch3idx; + std::vector muonCSCMatch1, muonCSCMatch1idx, muonCSCMatch2, muonCSCMatch2idx, muonCSCMatch3, muonCSCMatch3idx; + std::vector nSegments, nDTSegments, nCSCSegments; + + std::vector totPATmatches, totPATmatchesOS, totPATmatchesDisplID, totPATmatchesDisplIDOS, + totLoosePATmatchesDisplID, totLoosePATmatchesDisplIDOS; + std::vector LoosePATmatchesDisplIDOSpt, LoosePATmatchesDisplIDOSdsaDetID1, LoosePATmatchesDisplIDOSdsaDetID2, + LoosePATmatchesDisplIDOSdsaDetID3; + + for (unsigned int i = 0; i < nDSAMuons; i++) { + const reco::Track& dsaMuon = dsaMuons[i]; + idx.push_back(i); + + trkNumPlanes.push_back(dsaMuon.hitPattern().muonStationsWithValidHits()); + trkNumHits.push_back(dsaMuon.hitPattern().numberOfValidMuonHits()); + trkNumDTHits.push_back(dsaMuon.hitPattern().numberOfValidMuonDTHits()); + trkNumCSCHits.push_back(dsaMuon.hitPattern().numberOfValidMuonCSCHits()); + + outerEta.push_back(dsaMuon.extra().isNonnull() && dsaMuon.extra().isAvailable() ? dsaMuon.outerEta() : -999); + outerPhi.push_back(dsaMuon.extra().isNonnull() && dsaMuon.extra().isAvailable() ? dsaMuon.outerPhi() : -999); + + dzPV.push_back(dsaMuon.dz(pv.position())); + dzPVErr.push_back(std::hypot(dsaMuon.dzError(), pv.zError())); + reco::TransientTrack transientTrack = builder.build(dsaMuon); + TrajectoryStateClosestToPoint trajectoryPV = transientTrack.trajectoryStateClosestToPoint(primaryVertex); + dxyPVTraj.push_back(trajectoryPV.perigeeParameters().transverseImpactParameter()); + dxyPVTrajErr.push_back(trajectoryPV.perigeeError().transverseImpactParameterError()); + GlobalVector muonRefTrackDir(dsaMuon.px(), dsaMuon.py(), dsaMuon.pz()); + dxyPVSigned.push_back(IPTools::signedTransverseImpactParameter(transientTrack, muonRefTrackDir, pv).second.value()); + dxyPVSignedErr.push_back( + IPTools::signedTransverseImpactParameter(transientTrack, muonRefTrackDir, pv).second.error()); + + ip3DPVSigned.push_back(IPTools::signedImpactParameter3D(transientTrack, muonRefTrackDir, pv).second.value()); + ip3DPVSignedErr.push_back(IPTools::signedImpactParameter3D(transientTrack, muonRefTrackDir, pv).second.error()); + + float passesDisplacedId = 0; + if (passesDisplacedID(dsaMuon)) + passesDisplacedId = 1; + displacedId.push_back(passesDisplacedId); + + // Assigning 5 best matches and corresponding muon indices + std::vector> muonMatches(5, std::make_pair(-1.0, -1.0)); + std::vector> muonDTMatches(5, std::make_pair(-1.0, -1.0)); + std::vector> muonCSCMatches(5, std::make_pair(-1.0, -1.0)); + std::vector nMuonMatches; + std::vector nMuonDTMatches; + std::vector nMuonCSCMatches; + for (unsigned int j = 0; j < nMuons; j++) { + if (j > 4) + break; + const pat::Muon& muon = patMuons[j]; + // Muon-DSA Matches Table + auto [nDTMatches, nCSCMatches] = getMatches(muon, dsaMuon, minPositionDiffForMatching); + int nMatches = nDTMatches + nCSCMatches; + muonMatches[j] = std::make_pair(nMatches, j); + muonDTMatches[j] = std::make_pair(nDTMatches, j); + muonCSCMatches[j] = std::make_pair(nCSCMatches, j); + nMuonMatches.push_back(nMatches); + nMuonDTMatches.push_back(nDTMatches); + nMuonCSCMatches.push_back(nCSCMatches); + } + nMatchesPerMuon.push_back(nMuonMatches); + nMatchesPerMuon.push_back(nMuonDTMatches); + nMatchesPerMuon.push_back(nMuonCSCMatches); + std::sort(muonMatches.rbegin(), muonMatches.rend()); + std::sort(muonDTMatches.rbegin(), muonDTMatches.rend()); + std::sort(muonCSCMatches.rbegin(), muonCSCMatches.rend()); + muonMatch1.push_back(muonMatches[0].first); + muonMatch1idx.push_back(muonMatches[0].second); + muonMatch2.push_back(muonMatches[1].first); + muonMatch2idx.push_back(muonMatches[1].second); + muonMatch3.push_back(muonMatches[2].first); + muonMatch3idx.push_back(muonMatches[2].second); + muonMatch4.push_back(muonMatches[3].first); + muonMatch4idx.push_back(muonMatches[3].second); + muonMatch5.push_back(muonMatches[4].first); + muonMatch5idx.push_back(muonMatches[4].second); + + muonDTMatch1.push_back(muonDTMatches[0].first); + muonDTMatch1idx.push_back(muonDTMatches[0].second); + muonDTMatch2.push_back(muonDTMatches[1].first); + muonDTMatch2idx.push_back(muonDTMatches[1].second); + muonDTMatch3.push_back(muonDTMatches[2].first); + muonDTMatch3idx.push_back(muonDTMatches[2].second); + muonCSCMatch1.push_back(muonCSCMatches[0].first); + muonCSCMatch1idx.push_back(muonCSCMatches[0].second); + muonCSCMatch2.push_back(muonCSCMatches[1].first); + muonCSCMatch2idx.push_back(muonCSCMatches[1].second); + muonCSCMatch3.push_back(muonCSCMatches[2].first); + muonCSCMatch3idx.push_back(muonCSCMatches[2].second); + + auto [nDTSegments_, nCSCSegments_] = getSegments(dsaMuon); + nSegments.push_back(nDTSegments_ + nCSCSegments_); + nDTSegments.push_back(nDTSegments_); + nCSCSegments.push_back(nCSCSegments_); + } + + auto dsaMuonTab = std::make_unique(dsaMuons.size(), name_, false, true); + + dsaMuonTab->addColumn("idx", idx, ""); + dsaMuonTab->addColumn("trkNumPlanes", trkNumPlanes, ""); + dsaMuonTab->addColumn("trkNumHits", trkNumHits, ""); + dsaMuonTab->addColumn("trkNumDTHits", trkNumDTHits, ""); + dsaMuonTab->addColumn("trkNumCSCHits", trkNumCSCHits, ""); + + dsaMuonTab->addColumn("outerEta", outerEta, ""); + dsaMuonTab->addColumn("outerPhi", outerPhi, ""); + + dsaMuonTab->addColumn("dzPV", dzPV, ""); + dsaMuonTab->addColumn("dzPVErr", dzPVErr, ""); + dsaMuonTab->addColumn("dxyPVTraj", dxyPVTraj, ""); + dsaMuonTab->addColumn("dxyPVTrajErr", dxyPVTrajErr, ""); + dsaMuonTab->addColumn("dxyPVSigned", dxyPVSigned, ""); + dsaMuonTab->addColumn("dxyPVSignedErr", dxyPVSignedErr, ""); + dsaMuonTab->addColumn("ip3DPVSigned", ip3DPVSigned, ""); + dsaMuonTab->addColumn("ip3DPVSignedErr", ip3DPVSignedErr, ""); + + dsaMuonTab->addColumn("displacedID", displacedId, ""); + + dsaMuonTab->addColumn("muonMatch1", muonMatch1, ""); + dsaMuonTab->addColumn("muonMatch1idx", muonMatch1idx, ""); + dsaMuonTab->addColumn("muonMatch2", muonMatch2, ""); + dsaMuonTab->addColumn("muonMatch2idx", muonMatch2idx, ""); + dsaMuonTab->addColumn("muonMatch3", muonMatch3, ""); + dsaMuonTab->addColumn("muonMatch3idx", muonMatch3idx, ""); + dsaMuonTab->addColumn("muonMatch4", muonMatch4, ""); + dsaMuonTab->addColumn("muonMatch4idx", muonMatch4idx, ""); + dsaMuonTab->addColumn("muonMatch5", muonMatch5, ""); + dsaMuonTab->addColumn("muonMatch5idx", muonMatch5idx, ""); + + dsaMuonTab->addColumn("muonDTMatch1", muonDTMatch1, ""); + dsaMuonTab->addColumn("muonDTMatch1idx", muonDTMatch1idx, ""); + dsaMuonTab->addColumn("muonDTMatch2", muonDTMatch2, ""); + dsaMuonTab->addColumn("muonDTMatch2idx", muonDTMatch2idx, ""); + dsaMuonTab->addColumn("muonDTMatch3", muonDTMatch3, ""); + dsaMuonTab->addColumn("muonDTMatch3idx", muonDTMatch3idx, ""); + + dsaMuonTab->addColumn("muonCSCMatch1", muonCSCMatch1, ""); + dsaMuonTab->addColumn("muonCSCMatch1idx", muonCSCMatch1idx, ""); + dsaMuonTab->addColumn("muonCSCMatch2", muonCSCMatch2, ""); + dsaMuonTab->addColumn("muonCSCMatch2idx", muonCSCMatch2idx, ""); + dsaMuonTab->addColumn("muonCSCMatch3", muonCSCMatch3, ""); + dsaMuonTab->addColumn("muonCSCMatch3idx", muonCSCMatch3idx, ""); + + dsaMuonTab->addColumn("nSegments", nSegments, ""); + dsaMuonTab->addColumn("nDTSegments", nDTSegments, ""); + dsaMuonTab->addColumn("nCSCSegments", nCSCSegments, ""); + + iEvent.put(std::move(dsaMuonTab)); +} + +bool DSAMuonTableProducer::passesDisplacedID(const reco::Track& dsaMuon) const { + // displaced muon Id as recommended by Muon POG + float validHits = dsaMuon.hitPattern().numberOfValidMuonCSCHits() + dsaMuon.hitPattern().numberOfValidMuonDTHits(); + if (validHits > 12) { + if (dsaMuon.hitPattern().numberOfValidMuonCSCHits() != 0 || + (dsaMuon.hitPattern().numberOfValidMuonCSCHits() == 0 && dsaMuon.hitPattern().numberOfValidMuonDTHits() > 18)) { + if (dsaMuon.normalizedChi2() < 2.5) { + if (dsaMuon.ptError() / dsaMuon.pt() < 1) { + return true; + } + } + } + } + return false; +} + +// Returns number of DT and CSC segments of the DSA muon are associated to the PAT muon +std::tuple DSAMuonTableProducer::getMatches(const pat::Muon& muon, + const reco::Track& dsaMuon, + float minPositionDiff) const { + int nMatchesDT = 0; + int nMatchesCSC = 0; + + if (dsaMuon.extra().isNonnull() && dsaMuon.extra().isAvailable()) { + for (auto& hit : dsaMuon.recHits()) { + if (!hit->isValid()) + continue; + DetId id = hit->geographicalId(); + if (id.det() != DetId::Muon) + continue; + + if (id.subdetId() == MuonSubdetId::DT || id.subdetId() == MuonSubdetId::CSC) { + for (auto& chamber : muon.matches()) { + if (chamber.id.rawId() != id.rawId()) + continue; + + for (auto& segment : chamber.segmentMatches) { + if (fabs(segment.x - hit->localPosition().x()) < minPositionDiff && + fabs(segment.y - hit->localPosition().y()) < minPositionDiff) { + if (id.subdetId() == MuonSubdetId::DT) + nMatchesDT++; + else + nMatchesCSC++; + break; + } + } + } + } + } + } + return {nMatchesDT, nMatchesCSC}; +} + +std::tuple DSAMuonTableProducer::getSegments(const reco::Track& dsaMuon) const { + int nHitsDT = 0; + int nHitsCSC = 0; + if (dsaMuon.extra().isNonnull() && dsaMuon.extra().isAvailable()) { + for (auto& hit : dsaMuon.recHits()) { + if (!hit->isValid()) + continue; + DetId id = hit->geographicalId(); + if (id.det() != DetId::Muon) + continue; + if (id.subdetId() == MuonSubdetId::DT) + nHitsDT++; + if (id.subdetId() == MuonSubdetId::DT) + nHitsCSC++; + } + } + return {nHitsDT, nHitsCSC}; +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(DSAMuonTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc new file mode 100644 index 0000000000000..f8a27692b6b96 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc @@ -0,0 +1,338 @@ + +/* + +///////////////////////////////////// + +// Displaced tau training features with multi-threading : Mykyta Shchedrolosiev , Pritam Palit, created on 01/09/2025 + +///////////////////////////////////// + +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "FWCore/Utilities/interface/StreamID.h" +#include "PhysicsTools/TensorFlow/interface/TensorFlow.h" + +#include "PhysicsTools/NanoAOD/interface/DisTauTagScaling.h" + +void test_vector(std::vector& values) { + for (auto& value : values) { + if (std::isnan(value)) { + throw std::runtime_error("DisTauTag score output: NaN detected."); + } else if (std::isinf(value)) { + throw std::runtime_error("DisTauTag score output: Infinity detected."); + } else if (!std::isfinite(value)) { + throw std::runtime_error("DisTauTag score output: Non-standard value detected."); + } + } +} + +class DisTauTag : public edm::global::EDProducer<> { +public: + explicit DisTauTag(const edm::ParameterSet&); + ~DisTauTag() override; + + template + static Scalar getDeltaPhi(Scalar phi1, Scalar phi2); + + static void fill_zero(tensorflow::Tensor&); + +private: + void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override; + + template + float Scale(const Int_t, const Float_t, const bool) const; + void saveInputs(const tensorflow::Tensor& tensor, + const std::string& block_name, + const edm::EventID& eventId, + size_t jetIndex) const; + + void initializeTensorFlow(); + + const edm::FileInPath graphPath_; + const edm::EDGetTokenT jets_token; + const bool save_inputs_; + const unsigned batch_size_; + + // TensorFlow resources for threading + tensorflow::GraphDef* graphDef_; + tensorflow::Session* session_; +}; + +DisTauTag::~DisTauTag() { + if (session_) { + tensorflow::closeSession(session_); + } + delete graphDef_; +} + +DisTauTag::DisTauTag(const edm::ParameterSet& config) + : graphPath_(config.getParameter("graphPath")), + jets_token(consumes(config.getParameter("jets"))), + save_inputs_(config.getParameter("save_inputs")), + batch_size_(config.getParameter("batchSize")), + graphDef_(nullptr), + session_(nullptr) { + produces>("score0"); + produces>("score1"); + initializeTensorFlow(); +} + +void DisTauTag::initializeTensorFlow() { + if (!session_) { + const std::string graphFile = graphPath_.fullPath(); + graphDef_ = tensorflow::loadGraphDef(graphFile); + session_ = tensorflow::createSession(graphDef_); + } +} + +template +Scalar DisTauTag::getDeltaPhi(Scalar phi1, Scalar phi2) { + static constexpr Scalar pi = boost::math::constants::pi(); + Scalar dphi = phi1 - phi2; + if (dphi > pi) + dphi -= 2 * pi; + else if (dphi <= -pi) + dphi += 2 * pi; + return dphi; +} + +void DisTauTag::fill_zero(tensorflow::Tensor& tensor) { + size_t size_ = 1; + int num_dimensions = tensor.shape().dims(); + for (int ii_dim = 0; ii_dim < num_dimensions; ii_dim++) + size_ = size_ * tensor.shape().dim_size(ii_dim); + + for (size_t ii = 0; ii < size_; ii++) + tensor.flat()(ii) = 0.0; +} + +template +float DisTauTag::Scale(const Int_t idx, const Float_t value, const bool inner) const { + return std::clamp((value - FeatureT::mean.at(idx).at(inner)) / FeatureT::std.at(idx).at(inner), + FeatureT::lim_min.at(idx).at(inner), + FeatureT::lim_max.at(idx).at(inner)); +} + +void DisTauTag::saveInputs(const tensorflow::Tensor& tensor, + const std::string& block_name, + const edm::EventID& eventId, + size_t jetIndex) const { + const int tau_n = tensor.shape().dim_size(0); + const int pf_n = tensor.shape().dim_size(1); + const int ftr_n = tensor.shape().dim_size(2); + + std::string json_file_name = "distag_" + std::to_string(eventId.run()) + "_" + + std::to_string(eventId.luminosityBlock()) + "_" + std::to_string(eventId.event()) + "_" + + "jet_" + std::to_string(jetIndex) + ".json"; + + std::ofstream json_file(json_file_name.data()); + if (!json_file.is_open()) { + edm::LogError("DisTauTag") << "Failed to open file for saving inputs: " << json_file_name; + return; + } + + json_file << "{\"" << block_name << "\":["; + for (int tau_idx = 0; tau_idx < tau_n; tau_idx++) { + json_file << "["; + for (int pf_idx = 0; pf_idx < pf_n; pf_idx++) { + json_file << "["; + for (int ftr_idx = 0; ftr_idx < ftr_n; ftr_idx++) { + json_file << std::setprecision(7) << std::fixed << tensor.tensor()(tau_idx, pf_idx, ftr_idx); + if (ftr_idx < ftr_n - 1) + json_file << ", "; + } + json_file << "]"; + if (pf_idx < pf_n - 1) + json_file << ", "; + } + json_file << "]"; + if (tau_idx < tau_n - 1) + json_file << ", "; + } + json_file << "]}"; +} + +void DisTauTag::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& setup) const { + auto jets = event.getHandle(jets_token); + + const size_t jets_size = jets->size(); + + std::vector v_score0(jets_size, -9); + std::vector v_score1(jets_size, -9); + + // Processing jets in batches + const size_t n_batches = (jets_size + batch_size_ - 1) / batch_size_; + + for (size_t batch_idx = 0; batch_idx < n_batches; ++batch_idx) { + const size_t jet_start = batch_idx * batch_size_; + const size_t jet_end = std::min((batch_idx + 1) * batch_size_, jets_size); + const size_t current_batch_size = jet_end - jet_start; + + // batched inputs + tensorflow::Tensor input_1(tensorflow::DT_FLOAT, + {static_cast(current_batch_size), setupSizes::nSeq_PfCand, setupSizes::n_PfCand}); + tensorflow::Tensor input_2( + tensorflow::DT_FLOAT, + {static_cast(current_batch_size), setupSizes::nSeq_PfCand, setupSizes::n_PfCandCategorical}); + fill_zero(input_1); + fill_zero(input_2); + + // Filling the batches + for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) { + const size_t jetIndex = jet_start + batch_jet_idx; + const auto& jet = jets->at(jetIndex); + const auto& jet_p4 = jet.polarP4(); + + // Get jet daughters and sort by pt + const size_t nDaughters = jet.numberOfDaughters(); + std::vector indices(nDaughters); + std::iota(indices.begin(), indices.end(), 0); + std::sort(indices.begin(), indices.end(), [&](size_t a, size_t b) { + const auto& daughter_1 = jet.daughterPtr(a); + const auto& daughter_2 = jet.daughterPtr(b); + return daughter_1->polarP4().pt() > daughter_2->polarP4().pt(); + }); + + size_t daughter_idx = 0; + size_t tensor_idx = 0; + + while (tensor_idx < static_cast(setupSizes::nSeq_PfCand) && daughter_idx < nDaughters) { + const auto& jet_daughter = jet.daughterPtr(indices.at(daughter_idx)); + const auto daughter = dynamic_cast(jet_daughter.get()); + ++daughter_idx; + + if (!daughter) + continue; + + auto getVecRef = [&](tensorflow::Tensor& tensor, auto _fe, Float_t value) { + const int _feature_idx = static_cast(_fe); + if (_feature_idx < 0) + return; + tensor.tensor()(batch_jet_idx, tensor_idx, _feature_idx) = + Scale::scaler_type>(_feature_idx, value, false); + }; + + { // General features + typedef PfCand_Features Br; + getVecRef(input_1, Br::pfCand_valid, 1.0); + getVecRef(input_1, Br::pfCand_pt, static_cast(daughter->polarP4().pt())); + getVecRef(input_1, Br::pfCand_eta, static_cast(daughter->polarP4().eta())); + getVecRef(input_1, Br::pfCand_phi, static_cast(daughter->polarP4().phi())); + getVecRef(input_1, Br::pfCand_mass, static_cast(daughter->polarP4().mass())); + getVecRef(input_1, Br::pfCand_charge, static_cast(daughter->charge())); + getVecRef(input_1, Br::pfCand_puppiWeight, static_cast(daughter->puppiWeight())); + getVecRef(input_1, Br::pfCand_puppiWeightNoLep, static_cast(daughter->puppiWeightNoLep())); + getVecRef(input_1, Br::pfCand_lostInnerHits, static_cast(daughter->lostInnerHits())); + getVecRef(input_1, Br::pfCand_nPixelHits, static_cast(daughter->numberOfPixelHits())); + getVecRef(input_1, Br::pfCand_nHits, static_cast(daughter->numberOfHits())); + getVecRef(input_1, Br::pfCand_caloFraction, static_cast(daughter->caloFraction())); + getVecRef(input_1, Br::pfCand_hcalFraction, static_cast(daughter->hcalFraction())); + getVecRef(input_1, Br::pfCand_rawCaloFraction, static_cast(daughter->rawCaloFraction())); + getVecRef(input_1, Br::pfCand_rawHcalFraction, static_cast(daughter->rawHcalFraction())); + + getVecRef(input_1, Br::pfCand_hasTrackDetails, static_cast(daughter->hasTrackDetails())); + + if (daughter->hasTrackDetails()) { + if (std::isfinite(daughter->dz())) + getVecRef(input_1, Br::pfCand_dz, static_cast(daughter->dz())); + if (std::isfinite(daughter->dzError())) + getVecRef(input_1, Br::pfCand_dz_error, static_cast(daughter->dzError())); + if (std::isfinite(daughter->dxyError())) + getVecRef(input_1, Br::pfCand_dxy_error, static_cast(daughter->dxyError())); + + getVecRef(input_1, Br::pfCand_dxy, static_cast(daughter->dxy())); + getVecRef(input_1, Br::pfCand_track_chi2, static_cast(daughter->bestTrack()->chi2())); + getVecRef(input_1, Br::pfCand_track_ndof, static_cast(daughter->bestTrack()->ndof())); + } + + Float_t jet_eta = jet_p4.eta(); + Float_t jet_phi = jet_p4.phi(); + getVecRef(input_1, PfCand_Features::pfCand_deta, static_cast(daughter->polarP4().eta()) - jet_eta); + getVecRef(input_1, + PfCand_Features::pfCand_dphi, + getDeltaPhi(static_cast(daughter->polarP4().phi()), jet_phi)); + } + + { // Categorical features + typedef PfCandCategorical_Features Br; + getVecRef( + input_2, Br::pfCand_particleType, static_cast(TranslatePdgIdToPFParticleType(daughter->pdgId()))); + getVecRef(input_2, Br::pfCand_pvAssociationQuality, static_cast(daughter->pvAssociationQuality())); + getVecRef(input_2, Br::pfCand_fromPV, static_cast(daughter->fromPV())); + } + + ++tensor_idx; + } + } + + // Running inference on batch + std::vector outputs; + { tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); } + + // Storing results + for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) { + const size_t jetIndex = jet_start + batch_jet_idx; + v_score0[jetIndex] = outputs[0].matrix()(batch_jet_idx, 0); + v_score1[jetIndex] = outputs[0].matrix()(batch_jet_idx, 1); + } + + // Save inputs if requested (for debugging) + if (save_inputs_) { + for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) { + const size_t jetIndex = jet_start + batch_jet_idx; + + // Extract single jet inputs from the batch + tensorflow::Tensor single_input_1(tensorflow::DT_FLOAT, {1, setupSizes::nSeq_PfCand, setupSizes::n_PfCand}); + tensorflow::Tensor single_input_2(tensorflow::DT_FLOAT, + {1, setupSizes::nSeq_PfCand, setupSizes::n_PfCandCategorical}); + + for (size_t i = 0; i < static_cast(setupSizes::nSeq_PfCand); ++i) { + for (size_t j = 0; j < static_cast(setupSizes::n_PfCand); ++j) { + single_input_1.tensor()(0, i, j) = input_1.tensor()(batch_jet_idx, i, j); + } + for (size_t j = 0; j < static_cast(setupSizes::n_PfCandCategorical); ++j) { + single_input_2.tensor()(0, i, j) = input_2.tensor()(batch_jet_idx, i, j); + } + } + + saveInputs(single_input_1, "PfCand", event.id(), jetIndex); + saveInputs(single_input_2, "PfCandCategorical", event.id(), jetIndex); + } + } + } + + test_vector(v_score0); + test_vector(v_score1); + + std::unique_ptr> vm_score0(new edm::ValueMap()); + edm::ValueMap::Filler filler_score0(*vm_score0); + filler_score0.insert(jets, v_score0.begin(), v_score0.end()); + filler_score0.fill(); + event.put(std::move(vm_score0), "score0"); + + std::unique_ptr> vm_score1(new edm::ValueMap()); + edm::ValueMap::Filler filler_score1(*vm_score1); + filler_score1.insert(jets, v_score1.begin(), v_score1.end()); + filler_score1.fill(); + event.put(std::move(vm_score1), "score1"); +} + +DEFINE_FWK_MODULE(DisTauTag); diff --git a/PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc b/PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc new file mode 100644 index 0000000000000..07972b6548991 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc @@ -0,0 +1,491 @@ +// -*- C++ -*- +// +// Package: PhysicsTools/EXOnanoAOD +// Class: DispJetTableProducer +// +/**\class DispJetTableProducer + + Description: Additional variables for displaced vertices involving at least one lepton and other tracks + +*/ +// +// Original Author: Kirill Skovpen +// Created: Sat, 1 Mar 2025 08:37:01 GMT +// + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/MuonReco/interface/MuonSelectors.h" + +// nanoAOD include files +#include "DataFormats/NanoAOD/interface/FlatTable.h" + +// object specific include files +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/IPTools/interface/IPTools.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +class DispJetTableProducer : public edm::stream::EDProducer<> { +protected: + edm::ESGetToken ttbToken_; + edm::EDGetTokenT> electronTag_; + edm::EDGetTokenT> muonTag_; + edm::EDGetTokenT vtxTag_; + edm::EDGetTokenT secVtxTag_; + +public: + DispJetTableProducer(edm::ParameterSet const& params) + : ttbToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))), + electronTag_(consumes>(params.getParameter("electrons"))), + muonTag_(consumes>(params.getParameter("muons"))), + vtxTag_(consumes(params.getParameter("primaryVertex"))), + secVtxTag_(consumes(params.getParameter("secondaryVertex"))) { + produces("DispJetElectron"); + produces("DispJetElectronVtx"); + produces("DispJetElectronTrk"); + produces("DispJetMuon"); + produces("DispJetMuonVtx"); + produces("DispJetMuonTrk"); + } + + ~DispJetTableProducer() override {} + + void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override { + const TransientTrackBuilder* theB = &iSetup.getData(ttbToken_); + + const std::vector& electrons = iEvent.get(electronTag_); + const std::vector& muons = iEvent.get(muonTag_); + const reco::VertexCollection& primaryVertices = iEvent.get(vtxTag_); + const reco::VertexCollection& secVertices = iEvent.get(secVtxTag_); + const auto& pv = primaryVertices.at(0); + GlobalPoint primaryVertex(pv.x(), pv.y(), pv.z()); + + unsigned int nElectrons = electrons.size(); + unsigned int nMuons = muons.size(); + + std::vector el_idx; + std::vector el_lIVF_match; + std::vector el_IVF_df, el_IVF_ntracks, el_IVF_elid; + std::vector el_IVF_x, el_IVF_y, el_IVF_z, el_IVF_cx, el_IVF_cy, el_IVF_cz, el_IVF_chi2, el_IVF_pt, + el_IVF_eta, el_IVF_phi, el_IVF_E, el_IVF_mass; + std::vector el_IVF_trackcharge, el_IVF_trackelid, el_IVF_trackvtxid; + std::vector el_IVF_trackpt, el_IVF_tracketa, el_IVF_trackphi, el_IVF_trackE, el_IVF_trackdxy, el_IVF_trackdz; + std::vector el_IVF_tracksignedIP2D, el_IVF_tracksignedIP3D, el_IVF_tracksignedIP2Dsig, + el_IVF_tracksignedIP3Dsig; + std::vector el_IVF_signedIP2D, el_IVF_signedIP3D, el_IVF_signedIP2Dsig, el_IVF_signedIP3Dsig; + + std::vector mu_idx; + std::vector mu_lIVF_match; + std::vector mu_IVF_df, mu_IVF_ntracks, mu_IVF_muid; + std::vector mu_IVF_x, mu_IVF_y, mu_IVF_z, mu_IVF_cx, mu_IVF_cy, mu_IVF_cz, mu_IVF_chi2, mu_IVF_pt, + mu_IVF_eta, mu_IVF_phi, mu_IVF_E, mu_IVF_mass; + std::vector mu_IVF_trackcharge, mu_IVF_trackmuid, mu_IVF_trackvtxid; + std::vector mu_IVF_trackpt, mu_IVF_tracketa, mu_IVF_trackphi, mu_IVF_trackE, mu_IVF_trackdxy, mu_IVF_trackdz; + std::vector mu_IVF_tracksignedIP2D, mu_IVF_tracksignedIP3D, mu_IVF_tracksignedIP2Dsig, + mu_IVF_tracksignedIP3Dsig; + std::vector mu_IVF_signedIP2D, mu_IVF_signedIP3D, mu_IVF_signedIP2Dsig, mu_IVF_signedIP3Dsig; + + el_idx.clear(); + el_lIVF_match.clear(); + el_IVF_df.clear(); + el_IVF_ntracks.clear(); + el_IVF_elid.clear(); + el_IVF_x.clear(); + el_IVF_y.clear(); + el_IVF_z.clear(); + el_IVF_cx.clear(); + el_IVF_cy.clear(); + el_IVF_cz.clear(); + el_IVF_chi2.clear(); + el_IVF_pt.clear(); + el_IVF_eta.clear(); + el_IVF_phi.clear(); + el_IVF_E.clear(); + el_IVF_mass.clear(); + el_IVF_trackcharge.clear(); + el_IVF_trackelid.clear(); + el_IVF_trackvtxid.clear(); + el_IVF_trackpt.clear(); + el_IVF_tracketa.clear(); + el_IVF_trackphi.clear(); + el_IVF_trackE.clear(); + el_IVF_trackdxy.clear(); + el_IVF_trackdz.clear(); + el_IVF_tracksignedIP2D.clear(); + el_IVF_tracksignedIP2Dsig.clear(); + el_IVF_tracksignedIP3D.clear(); + el_IVF_tracksignedIP3Dsig.clear(); + + el_IVF_signedIP2D.clear(); + el_IVF_signedIP2Dsig.clear(); + el_IVF_signedIP3D.clear(); + el_IVF_signedIP3Dsig.clear(); + + mu_idx.clear(); + mu_lIVF_match.clear(); + mu_IVF_df.clear(); + mu_IVF_ntracks.clear(); + mu_IVF_muid.clear(); + mu_IVF_x.clear(); + mu_IVF_y.clear(); + mu_IVF_z.clear(); + mu_IVF_cx.clear(); + mu_IVF_cy.clear(); + mu_IVF_cz.clear(); + mu_IVF_chi2.clear(); + mu_IVF_pt.clear(); + mu_IVF_eta.clear(); + mu_IVF_phi.clear(); + mu_IVF_E.clear(); + mu_IVF_mass.clear(); + mu_IVF_trackcharge.clear(); + mu_IVF_trackmuid.clear(); + mu_IVF_trackvtxid.clear(); + mu_IVF_trackpt.clear(); + mu_IVF_tracketa.clear(); + mu_IVF_trackphi.clear(); + mu_IVF_trackE.clear(); + mu_IVF_trackdxy.clear(); + mu_IVF_trackdz.clear(); + mu_IVF_tracksignedIP2D.clear(); + mu_IVF_tracksignedIP2Dsig.clear(); + mu_IVF_tracksignedIP3D.clear(); + mu_IVF_tracksignedIP3Dsig.clear(); + + mu_IVF_signedIP2D.clear(); + mu_IVF_signedIP2Dsig.clear(); + mu_IVF_signedIP3D.clear(); + mu_IVF_signedIP3Dsig.clear(); + + int ntrack_max = 100; + int nElectronsSel = 0; + int nMuonsSel = 0; + + for (unsigned int i = 0; i < nElectrons; i++) { + const pat::Electron& el = electrons[i]; + + if (el.gsfTrack().isNull()) + continue; + if (el.pt() < 7) + continue; + if (fabs(el.eta()) > 2.5) + continue; + + el_idx.push_back(i); + el_lIVF_match.push_back(false); + + bool new_vtx = false; + double dR, deta, normchi2; + double mindR = 20, minnormchi2 = 10000; + int nVtx = 0; + reco::Vertex* vtxDisp = nullptr; + for (const reco::Vertex& vtx : secVertices) { + for (reco::Vertex::trackRef_iterator vtxTrackref = vtx.tracks_begin(); vtxTrackref != vtx.tracks_end(); + vtxTrackref++) { + reco::TrackRef vtxTrack = vtxTrackref->castTo(); + for (edm::Ref cand : el.associatedPackedPFCandidates()) { + dR = reco::deltaR(cand->eta(), cand->phi(), vtxTrack->eta(), vtxTrack->phi()); + deta = fabs(cand->eta() - vtxTrack->eta()); + normchi2 = fabs(vtx.chi2() / vtx.ndof()); + + if ((dR < 0.05 or (dR < 0.1 and deta < 0.03)) and + (dR < mindR or (dR == mindR and normchi2 < minnormchi2))) { + new_vtx = true; + vtxDisp = const_cast(&vtx); + el_lIVF_match[nElectronsSel] = true; + mindR = dR; + minnormchi2 = normchi2; + } + } + } + + if (new_vtx) { + el_IVF_x.push_back(vtx.x()); + el_IVF_y.push_back(vtx.y()); + el_IVF_z.push_back(vtx.z()); + el_IVF_cx.push_back(vtx.xError()); + el_IVF_cy.push_back(vtx.yError()); + el_IVF_cz.push_back(vtx.zError()); + el_IVF_df.push_back(vtx.ndof()); + el_IVF_chi2.push_back(vtx.chi2()); + el_IVF_pt.push_back(vtx.p4().pt()); + el_IVF_eta.push_back(vtx.p4().eta()); + el_IVF_phi.push_back(vtx.p4().phi()); + el_IVF_E.push_back(vtx.p4().energy()); + el_IVF_mass.push_back(vtx.p4().mass()); + el_IVF_elid.push_back(nElectronsSel); + + el_IVF_ntracks.push_back(0); + for (reco::Vertex::trackRef_iterator vtxTrackref = vtx.tracks_begin(); vtxTrackref != vtx.tracks_end(); + vtxTrackref++) { + if (el_IVF_ntracks.back() == ntrack_max) + break; + reco::TrackRef vtxTrack = vtxTrackref->castTo(); + const auto& trk = theB->build(vtxTrack); + Global3DVector dir(vtxTrack->px(), vtxTrack->py(), vtxTrack->pz()); + const auto& ip2d = IPTools::signedTransverseImpactParameter(trk, dir, *vtxDisp); + const auto& ip3d = IPTools::signedImpactParameter3D(trk, dir, *vtxDisp); + el_IVF_tracksignedIP2D.push_back(ip2d.second.value()); + el_IVF_tracksignedIP3D.push_back(ip3d.second.value()); + el_IVF_tracksignedIP2Dsig.push_back(ip2d.second.significance()); + el_IVF_tracksignedIP3Dsig.push_back(ip3d.second.significance()); + el_IVF_trackpt.push_back(vtxTrack->pt()); + el_IVF_tracketa.push_back(vtxTrack->eta()); + el_IVF_trackphi.push_back(vtxTrack->phi()); + el_IVF_trackE.push_back(vtxTrack->p()); + el_IVF_trackcharge.push_back(vtxTrack->charge()); + el_IVF_trackdxy.push_back(std::abs(vtxTrack->dxy(pv.position()))); + el_IVF_trackdz.push_back(std::abs(vtxTrack->dz(pv.position()))); + el_IVF_trackelid.push_back(nElectronsSel); + el_IVF_trackvtxid.push_back(nVtx); + el_IVF_ntracks.back()++; + } + nVtx++; + new_vtx = false; + + auto track = el.gsfTrack(); + if (track.isNonnull()) { + const auto& trk = theB->build(track); + Global3DVector dir(track->px(), track->py(), track->pz()); + const auto& ip2d = IPTools::signedTransverseImpactParameter(trk, dir, *vtxDisp); + const auto& ip3d = IPTools::signedImpactParameter3D(trk, dir, *vtxDisp); + el_IVF_signedIP2D.push_back(ip2d.second.value()); + el_IVF_signedIP3D.push_back(ip3d.second.value()); + el_IVF_signedIP2Dsig.push_back(ip2d.second.significance()); + el_IVF_signedIP3Dsig.push_back(ip3d.second.significance()); + } else { + el_IVF_signedIP2D.push_back(-999); + el_IVF_signedIP3D.push_back(-999); + el_IVF_signedIP2Dsig.push_back(-999); + el_IVF_signedIP3Dsig.push_back(-999); + } + } + } + nElectronsSel += 1; + } + + for (unsigned int i = 0; i < nMuons; i++) { + const pat::Muon& mu = muons[i]; + + if (mu.innerTrack().isNull()) + continue; + if (mu.pt() < 5) + continue; + if (fabs(mu.eta()) > 2.4) + continue; + if (!mu.isPFMuon()) + continue; + if (!(mu.isTrackerMuon() || mu.isGlobalMuon())) + continue; + + mu_idx.push_back(i); + mu_lIVF_match.push_back(false); + + bool new_vtx = false; + double ptdiff, normchi2; + double minptdiff = 10, minnormchi2 = 10000; + int nVtx = 0; + reco::Vertex* vtxDisp = nullptr; + for (const reco::Vertex& vtx : secVertices) { + for (reco::Vertex::trackRef_iterator vtxTrackref = vtx.tracks_begin(); vtxTrackref != vtx.tracks_end(); + vtxTrackref++) { + reco::TrackRef vtxTrack = vtxTrackref->castTo(); + for (size_t iCand = 0; iCand < mu.numberOfSourceCandidatePtrs(); ++iCand) { + if (!(mu.sourceCandidatePtr(iCand).isNonnull() and mu.sourceCandidatePtr(iCand).isAvailable())) + continue; + ptdiff = fabs(mu.sourceCandidatePtr(iCand)->pt() - vtxTrack->pt()); + normchi2 = fabs(vtx.chi2() / vtx.ndof()); + + if (ptdiff < 0.001 and (ptdiff < minptdiff or (ptdiff == minptdiff and normchi2 < minnormchi2))) { + new_vtx = true; + vtxDisp = const_cast(&vtx); + mu_lIVF_match[nMuonsSel] = true; + minptdiff = ptdiff; + minnormchi2 = normchi2; + } + } + } + if (new_vtx) { + mu_IVF_x.push_back(vtx.x()); + mu_IVF_y.push_back(vtx.y()); + mu_IVF_z.push_back(vtx.z()); + mu_IVF_cx.push_back(vtx.xError()); + mu_IVF_cy.push_back(vtx.yError()); + mu_IVF_cz.push_back(vtx.zError()); + mu_IVF_df.push_back(vtx.ndof()); + mu_IVF_chi2.push_back(vtx.chi2()); + mu_IVF_pt.push_back(vtx.p4().pt()); + mu_IVF_eta.push_back(vtx.p4().eta()); + mu_IVF_phi.push_back(vtx.p4().phi()); + mu_IVF_E.push_back(vtx.p4().energy()); + mu_IVF_mass.push_back(vtx.p4().mass()); + mu_IVF_muid.push_back(nMuonsSel); + + mu_IVF_ntracks.push_back(0); + for (reco::Vertex::trackRef_iterator vtxTrackref = vtx.tracks_begin(); vtxTrackref != vtx.tracks_end(); + vtxTrackref++) { + if (mu_IVF_ntracks.back() == ntrack_max) + break; + reco::TrackRef vtxTrack = vtxTrackref->castTo(); + const auto& trk = theB->build(vtxTrack); + Global3DVector dir(vtxTrack->px(), vtxTrack->py(), vtxTrack->pz()); + const auto& ip2d = IPTools::signedTransverseImpactParameter(trk, dir, *vtxDisp); + const auto& ip3d = IPTools::signedImpactParameter3D(trk, dir, *vtxDisp); + mu_IVF_tracksignedIP2D.push_back(ip2d.second.value()); + mu_IVF_tracksignedIP3D.push_back(ip3d.second.value()); + mu_IVF_tracksignedIP2Dsig.push_back(ip2d.second.significance()); + mu_IVF_tracksignedIP3Dsig.push_back(ip3d.second.significance()); + mu_IVF_trackpt.push_back(vtxTrack->pt()); + mu_IVF_tracketa.push_back(vtxTrack->eta()); + mu_IVF_trackphi.push_back(vtxTrack->phi()); + mu_IVF_trackE.push_back(vtxTrack->p()); + mu_IVF_trackcharge.push_back(vtxTrack->charge()); + mu_IVF_trackdxy.push_back(std::abs(vtxTrack->dxy(pv.position()))); + mu_IVF_trackdz.push_back(std::abs(vtxTrack->dz(pv.position()))); + mu_IVF_trackmuid.push_back(nMuonsSel); + mu_IVF_trackvtxid.push_back(nVtx); + mu_IVF_ntracks.back()++; + } + nVtx++; + new_vtx = false; + + auto track = mu.innerTrack(); + if (track.isNonnull()) { + const auto& trk = theB->build(track); + Global3DVector dir(track->px(), track->py(), track->pz()); + const auto& ip2d = IPTools::signedTransverseImpactParameter(trk, dir, *vtxDisp); + const auto& ip3d = IPTools::signedImpactParameter3D(trk, dir, *vtxDisp); + mu_IVF_signedIP2D.push_back(ip2d.second.value()); + mu_IVF_signedIP3D.push_back(ip3d.second.value()); + mu_IVF_signedIP2Dsig.push_back(ip2d.second.significance()); + mu_IVF_signedIP3Dsig.push_back(ip3d.second.significance()); + } else { + mu_IVF_signedIP2D.push_back(-999); + mu_IVF_signedIP3D.push_back(-999); + mu_IVF_signedIP2Dsig.push_back(-999); + mu_IVF_signedIP3Dsig.push_back(-999); + } + } + } + nMuonsSel += 1; + } + + auto dispJetElectronTab = std::make_unique(nElectronsSel, "DispJetElectron", false, false); + auto dispJetMuonTab = std::make_unique(nMuonsSel, "DispJetMuon", false, false); + + dispJetElectronTab->addColumn("idx", el_idx, ""); + dispJetElectronTab->addColumn("lIVF_match", el_lIVF_match, ""); + + auto dispJetElectronVtxTab = + std::make_unique(el_IVF_x.size(), "DispJetElectronVtx", false, false); + dispJetElectronVtxTab->addColumn("IVF_df", el_IVF_df, ""); + dispJetElectronVtxTab->addColumn("IVF_ntracks", el_IVF_ntracks, ""); + dispJetElectronVtxTab->addColumn("IVF_elid", el_IVF_elid, ""); + dispJetElectronVtxTab->addColumn("IVF_x", el_IVF_x, ""); + dispJetElectronVtxTab->addColumn("IVF_y", el_IVF_y, ""); + dispJetElectronVtxTab->addColumn("IVF_z", el_IVF_z, ""); + dispJetElectronVtxTab->addColumn("IVF_cx", el_IVF_cx, ""); + dispJetElectronVtxTab->addColumn("IVF_cy", el_IVF_cy, ""); + dispJetElectronVtxTab->addColumn("IVF_cz", el_IVF_cz, ""); + dispJetElectronVtxTab->addColumn("IVF_chi2", el_IVF_chi2, ""); + dispJetElectronVtxTab->addColumn("IVF_pt", el_IVF_pt, ""); + dispJetElectronVtxTab->addColumn("IVF_eta", el_IVF_eta, ""); + dispJetElectronVtxTab->addColumn("IVF_phi", el_IVF_phi, ""); + dispJetElectronVtxTab->addColumn("IVF_E", el_IVF_E, ""); + dispJetElectronVtxTab->addColumn("IVF_mass", el_IVF_mass, ""); + dispJetElectronVtxTab->addColumn("IVF_signedIP2D", el_IVF_signedIP2D, ""); + dispJetElectronVtxTab->addColumn("IVF_signedIP2Dsig", el_IVF_signedIP2Dsig, ""); + dispJetElectronVtxTab->addColumn("IVF_signedIP3D", el_IVF_signedIP3D, ""); + dispJetElectronVtxTab->addColumn("IVF_signedIP3Dsig", el_IVF_signedIP3Dsig, ""); + + int nTracksElectron = 0; + for (unsigned int iv = 0; iv < el_IVF_ntracks.size(); iv++) { + nTracksElectron += std::min(el_IVF_ntracks[iv], ntrack_max); + } + auto dispJetElectronTrkTab = + std::make_unique(nTracksElectron, "DispJetElectronTrk", false, false); + dispJetElectronTrkTab->addColumn("IVF_trackcharge", el_IVF_trackcharge, ""); + dispJetElectronTrkTab->addColumn("IVF_trackpt", el_IVF_trackpt, ""); + dispJetElectronTrkTab->addColumn("IVF_tracketa", el_IVF_tracketa, ""); + dispJetElectronTrkTab->addColumn("IVF_trackphi", el_IVF_trackphi, ""); + dispJetElectronTrkTab->addColumn("IVF_trackE", el_IVF_trackE, ""); + dispJetElectronTrkTab->addColumn("IVF_trackdxy", el_IVF_trackdxy, ""); + dispJetElectronTrkTab->addColumn("IVF_trackdz", el_IVF_trackdz, ""); + dispJetElectronTrkTab->addColumn("IVF_trackelid", el_IVF_trackelid, ""); + dispJetElectronTrkTab->addColumn("IVF_trackvtxid", el_IVF_trackvtxid, ""); + dispJetElectronTrkTab->addColumn("IVF_tracksignedIP2D", el_IVF_tracksignedIP2D, ""); + dispJetElectronTrkTab->addColumn("IVF_tracksignedIP2Dsig", el_IVF_tracksignedIP2Dsig, ""); + dispJetElectronTrkTab->addColumn("IVF_tracksignedIP3D", el_IVF_tracksignedIP3D, ""); + dispJetElectronTrkTab->addColumn("IVF_tracksignedIP3Dsig", el_IVF_tracksignedIP3Dsig, ""); + + dispJetMuonTab->addColumn("idx", mu_idx, ""); + dispJetMuonTab->addColumn("lIVF_match", mu_lIVF_match, ""); + + auto dispJetMuonVtxTab = std::make_unique(mu_IVF_x.size(), "DispJetMuonVtx", false, false); + dispJetMuonVtxTab->addColumn("IVF_df", mu_IVF_df, ""); + dispJetMuonVtxTab->addColumn("IVF_ntracks", mu_IVF_ntracks, ""); + dispJetMuonVtxTab->addColumn("IVF_muid", mu_IVF_muid, ""); + dispJetMuonVtxTab->addColumn("IVF_x", mu_IVF_x, ""); + dispJetMuonVtxTab->addColumn("IVF_y", mu_IVF_y, ""); + dispJetMuonVtxTab->addColumn("IVF_z", mu_IVF_z, ""); + dispJetMuonVtxTab->addColumn("IVF_cx", mu_IVF_cx, ""); + dispJetMuonVtxTab->addColumn("IVF_cy", mu_IVF_cy, ""); + dispJetMuonVtxTab->addColumn("IVF_cz", mu_IVF_cz, ""); + dispJetMuonVtxTab->addColumn("IVF_chi2", mu_IVF_chi2, ""); + dispJetMuonVtxTab->addColumn("IVF_pt", mu_IVF_pt, ""); + dispJetMuonVtxTab->addColumn("IVF_eta", mu_IVF_eta, ""); + dispJetMuonVtxTab->addColumn("IVF_phi", mu_IVF_phi, ""); + dispJetMuonVtxTab->addColumn("IVF_E", mu_IVF_E, ""); + dispJetMuonVtxTab->addColumn("IVF_mass", mu_IVF_mass, ""); + dispJetMuonVtxTab->addColumn("IVF_signedIP2D", mu_IVF_signedIP2D, ""); + dispJetMuonVtxTab->addColumn("IVF_signedIP2Dsig", mu_IVF_signedIP2Dsig, ""); + dispJetMuonVtxTab->addColumn("IVF_signedIP3D", mu_IVF_signedIP3D, ""); + dispJetMuonVtxTab->addColumn("IVF_signedIP3Dsig", mu_IVF_signedIP3Dsig, ""); + + int nTracksMuon = 0; + for (unsigned int iv = 0; iv < mu_IVF_ntracks.size(); iv++) { + nTracksMuon += std::min(mu_IVF_ntracks[iv], ntrack_max); + } + auto dispJetMuonTrkTab = std::make_unique(nTracksMuon, "DispJetMuonTrk", false, false); + dispJetMuonTrkTab->addColumn("IVF_trackcharge", mu_IVF_trackcharge, ""); + dispJetMuonTrkTab->addColumn("IVF_trackpt", mu_IVF_trackpt, ""); + dispJetMuonTrkTab->addColumn("IVF_tracketa", mu_IVF_tracketa, ""); + dispJetMuonTrkTab->addColumn("IVF_trackphi", mu_IVF_trackphi, ""); + dispJetMuonTrkTab->addColumn("IVF_trackE", mu_IVF_trackE, ""); + dispJetMuonTrkTab->addColumn("IVF_trackdxy", mu_IVF_trackdxy, ""); + dispJetMuonTrkTab->addColumn("IVF_trackdz", mu_IVF_trackdz, ""); + dispJetMuonTrkTab->addColumn("IVF_trackmuid", mu_IVF_trackmuid, ""); + dispJetMuonTrkTab->addColumn("IVF_trackvtxid", mu_IVF_trackvtxid, ""); + dispJetMuonTrkTab->addColumn("IVF_tracksignedIP2D", mu_IVF_tracksignedIP2D, ""); + dispJetMuonTrkTab->addColumn("IVF_tracksignedIP2Dsig", mu_IVF_tracksignedIP2Dsig, ""); + dispJetMuonTrkTab->addColumn("IVF_tracksignedIP3D", mu_IVF_tracksignedIP3D, ""); + dispJetMuonTrkTab->addColumn("IVF_tracksignedIP3Dsig", mu_IVF_tracksignedIP3Dsig, ""); + + iEvent.put(std::move(dispJetElectronTab), "DispJetElectron"); + iEvent.put(std::move(dispJetElectronVtxTab), "DispJetElectronVtx"); + iEvent.put(std::move(dispJetElectronTrkTab), "DispJetElectronTrk"); + iEvent.put(std::move(dispJetMuonTab), "DispJetMuon"); + iEvent.put(std::move(dispJetMuonVtxTab), "DispJetMuonVtx"); + iEvent.put(std::move(dispJetMuonTrkTab), "DispJetMuonTrk"); + } +}; + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(DispJetTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc b/PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc new file mode 100644 index 0000000000000..61ee491eb4a77 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc @@ -0,0 +1,295 @@ +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/NanoAOD/interface/FlatTable.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/IPTools/interface/IPTools.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +class ElectronExtendedTableProducer : public edm::global::EDProducer<> { +public: + explicit ElectronExtendedTableProducer(const edm::ParameterSet& iConfig) + : name_(iConfig.getParameter("name")), + rhoTag_(consumes(iConfig.getParameter("rho"))), + electronTag_(consumes>(iConfig.getParameter("electrons"))), + vtxTag_(consumes(iConfig.getParameter("primaryVertex"))), + jetTag_(consumes>(iConfig.getParameter("jets"))), + jetFatTag_(consumes>(iConfig.getParameter("jetsFat"))), + jetSubTag_(consumes>(iConfig.getParameter("jetsSub"))) { + produces(); + } + + ~ElectronExtendedTableProducer() override {}; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("rho")->setComment("input rho parameter"); + desc.add("electrons")->setComment("input electron collection"); + desc.add("primaryVertex")->setComment("input primary vertex collection"); + desc.add("jets")->setComment("input jet collection"); + desc.add("jetsFat")->setComment("input fat jet collection"); + desc.add("jetsSub")->setComment("input sub jet collection"); + desc.add("name")->setComment("name of the electron nanoaod::FlatTable we are extending"); + descriptions.add("electronTable", desc); + } + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + float GetDEtaInSeed(const pat::Electron* el) const; + float getPFIso(const pat::Electron& electron) const; + int findMatchedJet(const reco::Candidate& lepton, const std::vector& jets) const; + void fillLeptonJetVariables(const reco::GsfElectron* el, + const std::vector& jets, + const reco::Vertex& vertex, + const double rho, + std::vector* jetIdx, + std::vector relIso0p4, + std::vector* jetPtRatio, + std::vector* jetPtRel, + std::vector* jrtSelectedChargedMultiplicity) const; + + std::string name_; + edm::EDGetTokenT rhoTag_; + edm::EDGetTokenT> electronTag_; + edm::EDGetTokenT vtxTag_; + edm::EDGetTokenT> jetTag_; + edm::EDGetTokenT> jetFatTag_; + edm::EDGetTokenT> jetSubTag_; +}; + +void ElectronExtendedTableProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + const double rho = iEvent.get(rhoTag_); + const std::vector& electrons = iEvent.get(electronTag_); + const reco::VertexCollection& primaryVertices = iEvent.get(vtxTag_); + const std::vector& jets = iEvent.get(jetTag_); + const std::vector& jetFat = iEvent.get(jetFatTag_); + const std::vector& jetSub = iEvent.get(jetSubTag_); + + const auto& pv = primaryVertices.at(0); + GlobalPoint primaryVertex(pv.x(), pv.y(), pv.z()); + + unsigned int nElectrons = electrons.size(); + + std::vector idx, charge; + + std::vector relIso0p4; + std::vector jetPtRatio, jetPtRel; + std::vector jetIdx; + std::vector jetFatIdx, jetSubIdx; + std::vector jetSelectedChargedMultiplicity; + std::vector dxy, dz, IP3d, IP3dSig; + + std::vector isEB, isEE; + std::vector superClusterOverP, ecalEnergy, dEtaInSeed; + std::vector numberInnerHitsMissing, numberOfValidPixelHits, numberOfValidTrackerHits; + std::vector sigmaIetaIeta, deltaPhiSuperClusterTrack, deltaEtaSuperClusterTrack, eInvMinusPInv, hOverE; + + for (unsigned int i = 0; i < nElectrons; i++) { + const pat::Electron& electron = electrons[i]; + + idx.push_back(i); + + charge.push_back(electron.charge()); + + isEB.push_back(electron.isEB()); + isEE.push_back(electron.isEE()); + superClusterOverP.push_back(electron.eSuperClusterOverP()); + ecalEnergy.push_back(electron.ecalEnergy()); + dEtaInSeed.push_back(std::abs(GetDEtaInSeed(&electron))); + numberInnerHitsMissing.push_back( + electron.gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS)); + numberOfValidPixelHits.push_back( + (!electron.gsfTrack().isNull()) ? electron.gsfTrack()->hitPattern().numberOfValidPixelHits() : 0); + numberOfValidTrackerHits.push_back( + (!electron.gsfTrack().isNull()) ? electron.gsfTrack()->hitPattern().numberOfValidTrackerHits() : 0); + + sigmaIetaIeta.push_back(electron.full5x5_sigmaIetaIeta()); + deltaPhiSuperClusterTrack.push_back(fabs(electron.deltaPhiSuperClusterTrackAtVtx())); + deltaEtaSuperClusterTrack.push_back(fabs(electron.deltaEtaSuperClusterTrackAtVtx())); + eInvMinusPInv.push_back((1.0 - electron.eSuperClusterOverP()) / electron.correctedEcalEnergy()); + hOverE.push_back(electron.hadronicOverEm()); + + dxy.push_back(electron.dB(pat::Electron::PV2D)); + dz.push_back(electron.dB(pat::Electron::PVDZ)); + IP3d.push_back(electron.dB(pat::Electron::PV3D)); + IP3dSig.push_back(fabs(electron.dB(pat::Electron::PV3D) / electron.edB(pat::Electron::PV3D))); + + relIso0p4.push_back(getPFIso(electron)); + + fillLeptonJetVariables( + &electron, jets, pv, rho, &jetIdx, relIso0p4, &jetPtRatio, &jetPtRel, &jetSelectedChargedMultiplicity); + + const reco::Candidate* el_cand = dynamic_cast(&electron); + jetFatIdx.push_back(findMatchedJet(*el_cand, jetFat)); + jetSubIdx.push_back(findMatchedJet(*el_cand, jetSub)); + } + + auto tab = std::make_unique(nElectrons, name_, false, true); + tab->addColumn("idx", idx, "LLPnanoAOD electron index"); + + tab->addColumn("dxy", dxy, ""); + tab->addColumn("dz", dz, ""); + tab->addColumn("IP3d", IP3d, ""); + tab->addColumn("IP3dSig", IP3dSig, ""); + + tab->addColumn("relIso0p4", relIso0p4, ""); + tab->addColumn("jetPtRatio", jetPtRatio, ""); + tab->addColumn("jetPtRel", jetPtRel, ""); + tab->addColumn("jetSelectedChargedMultiplicity", jetSelectedChargedMultiplicity, ""); + tab->addColumn("jetIdx", jetIdx, ""); + tab->addColumn("jetFatIdx", jetFatIdx, ""); + tab->addColumn("jetSubIdx", jetSubIdx, ""); + + tab->addColumn("isEB", isEB, ""); + tab->addColumn("isEE", isEE, ""); + tab->addColumn("superClusterOverP", superClusterOverP, ""); + tab->addColumn("ecalEnergy", ecalEnergy, ""); + tab->addColumn("dEtaInSeed", dEtaInSeed, ""); + tab->addColumn("numberInnerHitsMissing", numberInnerHitsMissing, ""); + tab->addColumn("numberOfValidPixelHits", numberOfValidPixelHits, ""); + tab->addColumn("numberOfValidTrackerHits", numberOfValidTrackerHits, ""); + tab->addColumn("sigmaIetaIeta", sigmaIetaIeta, ""); + tab->addColumn("deltaPhiSuperClusterTrack", deltaPhiSuperClusterTrack, ""); + tab->addColumn("deltaEtaSuperClusterTrack", deltaEtaSuperClusterTrack, ""); + tab->addColumn("eInvMinusPInv", eInvMinusPInv, ""); + tab->addColumn("hOverE", hOverE, ""); + + iEvent.put(std::move(tab)); +} + +float ElectronExtendedTableProducer::GetDEtaInSeed(const pat::Electron* el) const { + if (el->superCluster().isNonnull() and el->superCluster()->seed().isNonnull()) + return el->deltaEtaSuperClusterTrackAtVtx() - el->superCluster()->eta() + el->superCluster()->seed()->eta(); + else + return std::numeric_limits::max(); +} + +template +bool isSourceCandidatePtrMatch(const T1& lhs, const T2& rhs) { + for (size_t lhsIndex = 0; lhsIndex < lhs.numberOfSourceCandidatePtrs(); ++lhsIndex) { + auto lhsSourcePtr = lhs.sourceCandidatePtr(lhsIndex); + for (size_t rhsIndex = 0; rhsIndex < rhs.numberOfSourceCandidatePtrs(); ++rhsIndex) { + auto rhsSourcePtr = rhs.sourceCandidatePtr(rhsIndex); + if (lhsSourcePtr == rhsSourcePtr) { + return true; + } + } + } + + return false; +} + +int ElectronExtendedTableProducer::findMatchedJet(const reco::Candidate& lepton, + const std::vector& jets) const { + int iJet = -1; + + unsigned int nJets = jets.size(); + + for (unsigned int i = 0; i < nJets; i++) { + const pat::Jet& jet = jets[i]; + if (isSourceCandidatePtrMatch(lepton, jet)) { + return i; + } + } + + return iJet; +} + +void ElectronExtendedTableProducer::fillLeptonJetVariables(const reco::GsfElectron* el, + const std::vector& jets, + const reco::Vertex& vertex, + const double rho, + std::vector* jetIdx, + std::vector relIso0p4, + std::vector* jetPtRatio, + std::vector* jetPtRel, + std::vector* jetSelectedChargedMultiplicity) const { + const reco::Candidate* cand = dynamic_cast(el); + int matchedJetIdx = findMatchedJet(*cand, jets); + + jetIdx->push_back(matchedJetIdx); + + if (matchedJetIdx < 0) { + float ptRatio = (1. / (1. + relIso0p4.back())); + jetPtRatio->push_back(ptRatio); + jetPtRel->push_back(0); + jetSelectedChargedMultiplicity->push_back(0); + } else { + const pat::Jet& jet = jets[matchedJetIdx]; + auto rawJetP4 = jet.correctedP4("Uncorrected"); + auto leptonP4 = cand->p4(); + + bool leptonEqualsJet = ((rawJetP4 - leptonP4).P() < 1e-4); + + if (leptonEqualsJet) { + jetPtRatio->push_back(1); + jetPtRel->push_back(0); + jetSelectedChargedMultiplicity->push_back(0); + } else { + auto L1JetP4 = jet.correctedP4("L1FastJet"); + double L2L3JEC = jet.pt() / L1JetP4.pt(); + auto lepAwareJetP4 = (L1JetP4 - leptonP4) * L2L3JEC + leptonP4; + + float ptRatio = cand->pt() / lepAwareJetP4.pt(); + float ptRel = leptonP4.Vect().Cross((lepAwareJetP4 - leptonP4).Vect().Unit()).R(); + jetPtRatio->push_back(ptRatio); + jetPtRel->push_back(ptRel); + jetSelectedChargedMultiplicity->push_back(0); + + for (const auto& daughterPtr : jet.daughterPtrVector()) { + const pat::PackedCandidate& daughter = *((const pat::PackedCandidate*)daughterPtr.get()); + + if (daughter.charge() == 0) + continue; + if (daughter.fromPV() < 2) + continue; + if (reco::deltaR(daughter, *cand) > 0.4) + continue; + if (!daughter.hasTrackDetails()) + continue; + + auto daughterTrack = daughter.pseudoTrack(); + + if (daughterTrack.pt() <= 1) + continue; + if (daughterTrack.hitPattern().numberOfValidHits() < 8) + continue; + if (daughterTrack.hitPattern().numberOfValidPixelHits() < 2) + continue; + if (daughterTrack.normalizedChi2() >= 5) + continue; + if (std::abs(daughterTrack.dz(vertex.position())) >= 17) + continue; + if (std::abs(daughterTrack.dxy(vertex.position())) >= 0.2) + continue; + ++jetSelectedChargedMultiplicity->back(); + } + } + } +} + +float ElectronExtendedTableProducer::getPFIso(const pat::Electron& electron) const { + return electron.userFloat("PFIsoAll04") / electron.pt(); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +//define this as a plug-in +DEFINE_FWK_MODULE(ElectronExtendedTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/ElectronVertexTableProducer.cc b/PhysicsTools/NanoAOD/plugins/ElectronVertexTableProducer.cc new file mode 100644 index 0000000000000..a194366fcf7d0 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/ElectronVertexTableProducer.cc @@ -0,0 +1,468 @@ +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/transform.h" +#include "DataFormats/NanoAOD/interface/FlatTable.h" + +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/GsfTrackReco/interface/GsfTrack.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "TLorentzVector.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/PatCandidates/interface/Electron.h" + +#include "DataFormats/Math/interface/deltaR.h" + +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" + +#include "RecoVertex/KinematicFitPrimitives/interface/ParticleMass.h" +#include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraint.h" +#include "RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h" +#include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h" +#include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h" +#include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h" +#include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h" +#include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h" + +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" + +#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h" +#include "PhysicsTools/RecoUtils/interface/CheckHitPattern.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" + +#include "TrackingTools/IPTools/interface/IPTools.h" + +#include +#include + +class ElectronVertexTableProducer : public edm::global::EDProducer<> { +public: + explicit ElectronVertexTableProducer(const edm::ParameterSet& iConfig) + : electronTag_(consumes>(iConfig.getParameter("electrons"))), + bsTag_(consumes(iConfig.getParameter("beamspot"))), + pvTag_(consumes(iConfig.getParameter("primaryVertex"))), + magneticFieldToken_(esConsumes()), + tkerGeomToken_(esConsumes()), + tkerTopoToken_(esConsumes()), + transientTrackBuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) { + produces("ElectronVertex"); + produces("ElectronVertexRefittedTracks"); + } + + ~ElectronVertexTableProducer() override {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("electrons")->setComment("input pat electrons collection"); + desc.add("beamspot")->setComment("input beamspot collection"); + desc.add("primaryVertex")->setComment("input primaryVertex collection"); + descriptions.add("electronVertexTables", desc); + } + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + std::pair getVxy(const reco::Vertex electronVertex) const; + std::pair getVxyz(const reco::Vertex electronVertex) const; + + std::tuple getDistanceBetweenElectronTracks( + const reco::Track& track1, const reco::Track& track2, const edm::ESHandle& magneticField) const; + + const edm::EDGetTokenT> electronTag_; + const edm::EDGetTokenT bsTag_; + const edm::EDGetTokenT pvTag_; + const edm::ESGetToken magneticFieldToken_; + const edm::ESGetToken tkerGeomToken_; + const edm::ESGetToken tkerTopoToken_; + const edm::ESGetToken transientTrackBuilderToken_; +}; + +void ElectronVertexTableProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + const std::vector& electrons = iEvent.get(electronTag_); + + const reco::BeamSpot& beamSpotInput = iEvent.get(bsTag_); + const auto& bs = beamSpotInput.position(); + GlobalPoint beamSpot(bs.x(), bs.y(), bs.z()); + reco::Vertex beamSpotVertex(beamSpotInput.position(), beamSpotInput.covariance3D()); + + const reco::VertexCollection& primaryVertices = iEvent.get(pvTag_); + const auto& pv = primaryVertices.at(0); + GlobalPoint primaryVertex(pv.x(), pv.y(), pv.z()); + + auto const& magneticField = &iSetup.getData(magneticFieldToken_); + + auto const& tkerGeom = &iSetup.getData(tkerGeomToken_); + auto const& tkerTopo = &iSetup.getData(tkerTopoToken_); + const TransientTrackBuilder& builder = iSetup.getData(transientTrackBuilderToken_); + + KalmanVertexFitter vertexFitter(true, true); + + int nElectronVertices = 0; + int refittedTrackIdx_counter = 0; + + std::vector vertexIsValid; + std::vector vxy, vxyz, vx, vy, vz, t, vxySigma, vxyzSigma, vxErr, vyErr, vzErr, tErr; + std::vector chi2, ndof, normChi2, dR, originalElectronIdx1, originalElectronIdx2, refittedTrackIdx1, + refittedTrackIdx2; + std::vector DCA, DCAstatus, DCAx, DCAy, DCAz; + std::vector hitsInFrontOfVert1, missHitsAfterVert1, hitsInFrontOfVert2, missHitsAfterVert2; + + std::vector refittedTrackIdx, refittedTrackOriginalIdx, refittedTrackPt, refittedTrackPtErr, refittedTrackPx, + refittedTrackPy, refittedTrackPz; + std::vector refittedTrackEta, refittedTrackEtaErr, refittedTrackPhi, refittedTrackPhiErr, refittedTrackCharge, + refittedTrackNormChi2, refittedTrackNdof, refittedTrackChi2; + std::vector refittedTrackDzPV, refittedTrackDzPVErr, refittedTrackDxyPVTraj, refittedTrackDxyPVTrajErr, + refittedTrackDxyPVSigned, refittedTrackDxyPVSignedErr, refittedTrackIp3DPVSigned, refittedTrackIp3DPVSignedErr; + std::vector refittedTrackDxyBS, refittedTrackDxyBSErr, refittedTrackDzBS, refittedTrackDzBSErr, + refittedTrackDxyBSTraj, refittedTrackDxyBSTrajErr, refittedTrackDxyBSSigned, refittedTrackDxyBSSignedErr, + refittedTrackIp3DBSSigned, refittedTrackIp3DBSSignedErr; + std::vector refittedVertMass, refittedVertPt, refittedVertEta, refittedVertPhi; + // pat electrons + for (size_t i = 0; i < electrons.size(); i++) { + const pat::Electron& electron_i = electrons.at(i); + + reco::GsfTrackRef trackRef_i = electron_i.gsfTrack(); + + const auto& electronTrack_i = trackRef_i.get(); + reco::TransientTrack electronTransientTrack_i = builder.build(electronTrack_i); + + // pat-pat electron vertex + for (size_t j = i + 1; j < electrons.size(); j++) { + const pat::Electron& electron_j = electrons.at(j); + + reco::GsfTrackRef trackRef_j = electron_j.gsfTrack(); + + const auto& electronTrack_j = trackRef_j.get(); + reco::TransientTrack electronTransientTrack_j = builder.build(electronTrack_j); + + std::vector electronTransientTracks{}; + electronTransientTracks.push_back(electronTransientTrack_i); + electronTransientTracks.push_back(electronTransientTrack_j); + + TransientVertex transientElectronVertex = vertexFitter.vertex(electronTransientTracks); + reco::Vertex electronVertex = reco::Vertex(transientElectronVertex); + + if (!transientElectronVertex.isValid()) + continue; + std::tuple distanceTuple = + getDistanceBetweenElectronTracks(*electronTrack_i, *electronTrack_j, magneticField); + // if dca status is good but dca is more than 15 cm + if (std::get<1>(distanceTuple) && std::get<0>(distanceTuple) > 15) + continue; + + nElectronVertices++; + + vertexIsValid.push_back(transientElectronVertex.isValid()); + std::pair vxy_ = getVxy(electronVertex); + vxy.push_back(vxy_.first); + vxySigma.push_back(vxy_.second); + std::pair vxyz_ = getVxyz(electronVertex); + vxyz.push_back(vxyz_.first); + vxyzSigma.push_back(vxyz_.second); + ndof.push_back(electronVertex.ndof()); + normChi2.push_back(electronVertex.normalizedChi2()); + chi2.push_back(electronVertex.chi2()); + vx.push_back(electronVertex.x()); + vy.push_back(electronVertex.y()); + vz.push_back(electronVertex.z()); + t.push_back(electronVertex.t()); + vxErr.push_back(electronVertex.xError()); + vyErr.push_back(electronVertex.yError()); + vzErr.push_back(electronVertex.zError()); + tErr.push_back(electronVertex.tError()); + + dR.push_back(reco::deltaR(electron_i, electron_j)); + originalElectronIdx1.push_back(i); + originalElectronIdx2.push_back(j); + + DCA.push_back(std::get<0>(distanceTuple)); + DCAstatus.push_back(std::get<1>(distanceTuple)); + DCAx.push_back(std::get<2>(distanceTuple).x()); + DCAy.push_back(std::get<2>(distanceTuple).y()); + DCAz.push_back(std::get<2>(distanceTuple).z()); + + CheckHitPattern checkHitPattern; + checkHitPattern.init(tkerTopo, *tkerGeom, builder); + CheckHitPattern::Result hitPattern_i = checkHitPattern(*electronTrack_i, transientElectronVertex.vertexState()); + hitsInFrontOfVert1.push_back(hitPattern_i.hitsInFrontOfVert); + missHitsAfterVert1.push_back(hitPattern_i.missHitsAfterVert); + + CheckHitPattern::Result hitPattern_j = checkHitPattern(*electronTrack_j, transientElectronVertex.vertexState()); + hitsInFrontOfVert2.push_back(hitPattern_j.hitsInFrontOfVert); + missHitsAfterVert2.push_back(hitPattern_j.missHitsAfterVert); + + reco::TransientTrack refittedTrack_i = transientElectronVertex.refittedTrack(electronTransientTrack_i); + reco::TransientTrack refittedTrack_j = transientElectronVertex.refittedTrack(electronTransientTrack_j); + refittedTrackOriginalIdx.push_back(i); + refittedTrackPt.push_back(refittedTrack_i.track().pt()); + refittedTrackPtErr.push_back(refittedTrack_i.track().ptError()); + refittedTrackPx.push_back(refittedTrack_i.track().px()); + refittedTrackPy.push_back(refittedTrack_i.track().py()); + refittedTrackPz.push_back(refittedTrack_i.track().pz()); + refittedTrackEta.push_back(refittedTrack_i.track().eta()); + refittedTrackEtaErr.push_back(refittedTrack_i.track().etaError()); + refittedTrackPhi.push_back(refittedTrack_i.track().phi()); + refittedTrackPhiErr.push_back(refittedTrack_i.track().phiError()); + refittedTrackCharge.push_back(refittedTrack_i.track().charge()); + refittedTrackNormChi2.push_back(refittedTrack_i.normalizedChi2()); + refittedTrackNdof.push_back(refittedTrack_i.ndof()); + refittedTrackChi2.push_back(refittedTrack_i.chi2()); + + refittedTrackDzPV.push_back(refittedTrack_i.track().dz(pv.position())); + refittedTrackDzPVErr.push_back(std::hypot(refittedTrack_i.track().dzError(), pv.zError())); + TrajectoryStateClosestToPoint trajectoryPV_i = refittedTrack_i.trajectoryStateClosestToPoint(primaryVertex); + refittedTrackDxyPVTraj.push_back(trajectoryPV_i.perigeeParameters().transverseImpactParameter()); + refittedTrackDxyPVTrajErr.push_back(trajectoryPV_i.perigeeError().transverseImpactParameterError()); + GlobalVector electronRefTrackDir_i( + refittedTrack_i.track().px(), refittedTrack_i.track().py(), refittedTrack_i.track().pz()); + refittedTrackDxyPVSigned.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_i, electronRefTrackDir_i, pv).second.value()); + refittedTrackDxyPVSignedErr.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_i, electronRefTrackDir_i, pv).second.error()); + refittedTrackIp3DPVSigned.push_back( + IPTools::signedImpactParameter3D(refittedTrack_i, electronRefTrackDir_i, pv).second.value()); + refittedTrackIp3DPVSignedErr.push_back( + IPTools::signedImpactParameter3D(refittedTrack_i, electronRefTrackDir_i, pv).second.error()); + refittedTrackDxyBS.push_back(refittedTrack_i.track().dxy(bs)); + refittedTrackDxyBSErr.push_back(std::hypot(refittedTrack_i.track().dxyError(), beamSpotVertex.zError())); + refittedTrackDzBS.push_back(refittedTrack_i.track().dz(bs)); + refittedTrackDzBSErr.push_back(std::hypot(refittedTrack_i.track().dzError(), beamSpotVertex.zError())); + TrajectoryStateClosestToBeamLine trajectoryBS_i = refittedTrack_i.stateAtBeamLine(); + refittedTrackDxyBSTraj.push_back(trajectoryBS_i.transverseImpactParameter().value()); + refittedTrackDxyBSTrajErr.push_back(trajectoryBS_i.transverseImpactParameter().error()); + refittedTrackDxyBSSigned.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_i, electronRefTrackDir_i, beamSpotVertex) + .second.value()); + refittedTrackDxyBSSignedErr.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_i, electronRefTrackDir_i, beamSpotVertex) + .second.error()); + refittedTrackIp3DBSSigned.push_back( + IPTools::signedImpactParameter3D(refittedTrack_i, electronRefTrackDir_i, beamSpotVertex).second.value()); + refittedTrackIp3DBSSignedErr.push_back( + IPTools::signedImpactParameter3D(refittedTrack_i, electronRefTrackDir_i, beamSpotVertex).second.error()); + refittedTrackIdx.push_back(refittedTrackIdx_counter); + refittedTrackIdx1.push_back(refittedTrackIdx_counter); + refittedTrackIdx_counter++; + + refittedTrackOriginalIdx.push_back(j); + refittedTrackPt.push_back(refittedTrack_j.track().pt()); + refittedTrackPtErr.push_back(refittedTrack_j.track().ptError()); + refittedTrackPx.push_back(refittedTrack_j.track().px()); + refittedTrackPy.push_back(refittedTrack_j.track().py()); + refittedTrackPz.push_back(refittedTrack_j.track().pz()); + refittedTrackEta.push_back(refittedTrack_j.track().eta()); + refittedTrackEtaErr.push_back(refittedTrack_j.track().etaError()); + refittedTrackPhi.push_back(refittedTrack_j.track().phi()); + refittedTrackPhiErr.push_back(refittedTrack_j.track().phiError()); + refittedTrackCharge.push_back(refittedTrack_j.track().charge()); + refittedTrackNormChi2.push_back(refittedTrack_j.normalizedChi2()); + refittedTrackNdof.push_back(refittedTrack_j.ndof()); + refittedTrackChi2.push_back(refittedTrack_j.chi2()); + + refittedTrackDzPV.push_back(refittedTrack_j.track().dz(pv.position())); + refittedTrackDzPVErr.push_back(std::hypot(refittedTrack_j.track().dzError(), pv.zError())); + TrajectoryStateClosestToPoint trajectoryPV_j = refittedTrack_j.trajectoryStateClosestToPoint(primaryVertex); + refittedTrackDxyPVTraj.push_back(trajectoryPV_j.perigeeParameters().transverseImpactParameter()); + refittedTrackDxyPVTrajErr.push_back(trajectoryPV_j.perigeeError().transverseImpactParameterError()); + GlobalVector electronRefTrackDir_j( + refittedTrack_j.track().px(), refittedTrack_j.track().py(), refittedTrack_j.track().pz()); + refittedTrackDxyPVSigned.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_j, electronRefTrackDir_j, pv).second.value()); + refittedTrackDxyPVSignedErr.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_j, electronRefTrackDir_j, pv).second.error()); + refittedTrackIp3DPVSigned.push_back( + IPTools::signedImpactParameter3D(refittedTrack_j, electronRefTrackDir_j, pv).second.value()); + refittedTrackIp3DPVSignedErr.push_back( + IPTools::signedImpactParameter3D(refittedTrack_j, electronRefTrackDir_j, pv).second.error()); + refittedTrackDxyBS.push_back(refittedTrack_j.track().dxy(bs)); + refittedTrackDxyBSErr.push_back(std::hypot(refittedTrack_j.track().dxyError(), beamSpotVertex.zError())); + refittedTrackDzBS.push_back(refittedTrack_j.track().dz(bs)); + refittedTrackDzBSErr.push_back(std::hypot(refittedTrack_j.track().dzError(), beamSpotVertex.zError())); + TrajectoryStateClosestToBeamLine trajectoryBS_j = refittedTrack_j.stateAtBeamLine(); + refittedTrackDxyBSTraj.push_back(trajectoryBS_j.transverseImpactParameter().value()); + refittedTrackDxyBSTrajErr.push_back(trajectoryBS_j.transverseImpactParameter().error()); + refittedTrackDxyBSSigned.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_j, electronRefTrackDir_j, beamSpotVertex) + .second.value()); + refittedTrackDxyBSSignedErr.push_back( + IPTools::signedTransverseImpactParameter(refittedTrack_j, electronRefTrackDir_j, beamSpotVertex) + .second.error()); + refittedTrackIp3DBSSigned.push_back( + IPTools::signedImpactParameter3D(refittedTrack_j, electronRefTrackDir_j, beamSpotVertex).second.value()); + refittedTrackIp3DBSSignedErr.push_back( + IPTools::signedImpactParameter3D(refittedTrack_j, electronRefTrackDir_j, beamSpotVertex).second.error()); + refittedTrackIdx.push_back(refittedTrackIdx_counter); + refittedTrackIdx2.push_back(refittedTrackIdx_counter); + refittedTrackIdx_counter++; + + // Perform kinematic fit to re-compute dielectron four-momentum (especially for better mass resolution) + KinematicParticleFactoryFromTransientTrack pFactory; + ParticleMass e_mass = 0.000511; + float e_sigma = 0.00000001; + float chi = 0.; + float ndf = 0.; + std::vector eleParticles; + eleParticles.push_back(pFactory.particle(electronTransientTracks[0], e_mass, chi, ndf, e_sigma)); + eleParticles.push_back(pFactory.particle(electronTransientTracks[1], e_mass, chi, ndf, e_sigma)); + KinematicParticleVertexFitter fitter; + try { + RefCountedKinematicTree vertexFitTree = fitter.fit(eleParticles); + if (vertexFitTree->isValid()) { + vertexFitTree->movePointerToTheTop(); + auto diele_part = vertexFitTree->currentParticle(); + auto diele_state = diele_part->currentState(); + auto daughters = vertexFitTree->daughterParticles(); + refittedVertMass.push_back(diele_state.mass()); + refittedVertPt.push_back(diele_state.globalMomentum().transverse()); + refittedVertEta.push_back(diele_state.globalMomentum().eta()); + refittedVertPhi.push_back(diele_state.globalMomentum().phi()); + } else { + refittedVertMass.push_back(-999.); + refittedVertPt.push_back(-999.); + refittedVertEta.push_back(-999.); + refittedVertPhi.push_back(-999.); + } + } catch (std::exception& ex) { + std::cout << "kinematic vertex fit failed!" << std::endl; + refittedVertMass.push_back(-999.); + refittedVertPt.push_back(-999.); + refittedVertEta.push_back(-999.); + refittedVertPhi.push_back(-999.); + } + } + } + + auto vertexTab = std::make_unique(nElectronVertices, "ElectronVertex", false, false); + auto refittedTracksTab = + std::make_unique(nElectronVertices * 2, "ElectronVertexRefittedTracks", false, false); + + vertexTab->addColumn("isValid", vertexIsValid, ""); + vertexTab->addColumn("vxy", vxy, ""); + vertexTab->addColumn("vxySigma", vxySigma, ""); + vertexTab->addColumn("vxyz", vxyz, ""); + vertexTab->addColumn("vxyzSigma", vxyzSigma, ""); + vertexTab->addColumn("chi2", chi2, ""); + vertexTab->addColumn("ndof", ndof, ""); + vertexTab->addColumn("normChi2", normChi2, ""); + vertexTab->addColumn("vx", vx, ""); + vertexTab->addColumn("vy", vy, ""); + vertexTab->addColumn("vz", vz, ""); + vertexTab->addColumn("t", t, ""); + vertexTab->addColumn("vxErr", vxErr, ""); + vertexTab->addColumn("vyErr", vyErr, ""); + vertexTab->addColumn("vzErr", vzErr, ""); + vertexTab->addColumn("tErr", tErr, ""); + vertexTab->addColumn("dR", dR, ""); + vertexTab->addColumn("originalElectronIdx1", originalElectronIdx1, ""); + vertexTab->addColumn("originalElectronIdx2", originalElectronIdx2, ""); + vertexTab->addColumn("dca", DCA, ""); + vertexTab->addColumn("dcaStatus", DCAstatus, ""); + vertexTab->addColumn("dcax", DCAx, ""); + vertexTab->addColumn("dcay", DCAy, ""); + vertexTab->addColumn("dcaz", DCAz, ""); + vertexTab->addColumn("hitsInFrontOfVert1", hitsInFrontOfVert1, ""); + vertexTab->addColumn("hitsInFrontOfVert2", hitsInFrontOfVert2, ""); + vertexTab->addColumn("missHitsAfterVert1", missHitsAfterVert1, ""); + vertexTab->addColumn("missHitsAfterVert2", missHitsAfterVert2, ""); + vertexTab->addColumn("refittedTrackIdx1", refittedTrackIdx1, ""); + vertexTab->addColumn("refittedTrackIdx2", refittedTrackIdx2, ""); + vertexTab->addColumn("massRefit", refittedVertMass, ""); + vertexTab->addColumn("ptRefit", refittedVertPt, ""); + vertexTab->addColumn("etaRefit", refittedVertEta, ""); + vertexTab->addColumn("phiRefit", refittedVertPhi, ""); + + iEvent.put(std::move(vertexTab), "ElectronVertex"); + + refittedTracksTab->addColumn("idx", refittedTrackIdx, ""); + refittedTracksTab->addColumn("originalElectronIdx", refittedTrackOriginalIdx, ""); + refittedTracksTab->addColumn("pt", refittedTrackPt, ""); + refittedTracksTab->addColumn("ptErr", refittedTrackPtErr, ""); + refittedTracksTab->addColumn("px", refittedTrackPx, ""); + refittedTracksTab->addColumn("py", refittedTrackPy, ""); + refittedTracksTab->addColumn("pz", refittedTrackPz, ""); + refittedTracksTab->addColumn("eta", refittedTrackEta, ""); + refittedTracksTab->addColumn("etaErr", refittedTrackEtaErr, ""); + refittedTracksTab->addColumn("phi", refittedTrackPhi, ""); + refittedTracksTab->addColumn("phiErr", refittedTrackPhiErr, ""); + refittedTracksTab->addColumn("charge", refittedTrackCharge, ""); + refittedTracksTab->addColumn("normChi2", refittedTrackNormChi2, ""); + refittedTracksTab->addColumn("ndof", refittedTrackNdof, ""); + refittedTracksTab->addColumn("chi2", refittedTrackChi2, ""); + refittedTracksTab->addColumn("dzPV", refittedTrackDzPV, ""); + refittedTracksTab->addColumn("dzPVErr", refittedTrackDzPVErr, ""); + refittedTracksTab->addColumn("dxyPVTraj", refittedTrackDxyPVTraj, ""); + refittedTracksTab->addColumn("dxyPVTrajErr", refittedTrackDxyPVTrajErr, ""); + refittedTracksTab->addColumn("dxyPVSigned", refittedTrackDxyPVSigned, ""); + refittedTracksTab->addColumn("dxyPVSignedErr", refittedTrackDxyPVSignedErr, ""); + refittedTracksTab->addColumn("ip3DPVSigned", refittedTrackIp3DPVSigned, ""); + refittedTracksTab->addColumn("ip3DPVSignedErr", refittedTrackIp3DPVSignedErr, ""); + refittedTracksTab->addColumn("dxyBS", refittedTrackDxyBS, ""); + refittedTracksTab->addColumn("dxyBSErr", refittedTrackDxyBSErr, ""); + refittedTracksTab->addColumn("dzBS", refittedTrackDzBS, ""); + refittedTracksTab->addColumn("dzBSErr", refittedTrackDzBSErr, ""); + refittedTracksTab->addColumn("dxyBSTraj", refittedTrackDxyBSTraj, ""); + refittedTracksTab->addColumn("dxyBSTrajErr", refittedTrackDxyBSTrajErr, ""); + refittedTracksTab->addColumn("dxyBSSigned", refittedTrackDxyBSSigned, ""); + refittedTracksTab->addColumn("dxyBSSignedErr", refittedTrackDxyBSSignedErr, ""); + refittedTracksTab->addColumn("ip3DBSSigned", refittedTrackIp3DBSSigned, ""); + refittedTracksTab->addColumn("ip3DBSSignedErr", refittedTrackIp3DBSSignedErr, ""); + + iEvent.put(std::move(refittedTracksTab), "ElectronVertexRefittedTracks"); +} + +std::pair ElectronVertexTableProducer::getVxy(const reco::Vertex electronVertex) const { + float vxy = sqrt(electronVertex.x() * electronVertex.x() + electronVertex.y() * electronVertex.y()); + float vxySigma = + (1 / vxy) * sqrt(electronVertex.x() * electronVertex.x() * electronVertex.xError() * electronVertex.xError() + + electronVertex.y() * electronVertex.y() * electronVertex.yError() * electronVertex.yError()); + return std::make_pair(vxy, vxySigma); +} + +std::pair ElectronVertexTableProducer::getVxyz(const reco::Vertex electronVertex) const { + float vxyz = sqrt(electronVertex.x() * electronVertex.x() + electronVertex.y() * electronVertex.y() + + electronVertex.z() * electronVertex.z()); + float vxyzSigma = + (1 / vxyz) * sqrt(electronVertex.x() * electronVertex.x() * electronVertex.xError() * electronVertex.xError() + + electronVertex.y() * electronVertex.y() * electronVertex.yError() * electronVertex.yError() + + electronVertex.z() * electronVertex.z() * electronVertex.zError() * electronVertex.zError()); + return std::make_pair(vxyz, vxyzSigma); +} + +/** +* Proximity between the electrons based on EXO-23-010 +* Getting Distance of Closest Approach between electron tracks using TwoTrackMinimumDistance +* Returns tuple of distance (float), error of distance (float) and crossing point (GlobalPoint) +**/ +std::tuple ElectronVertexTableProducer::getDistanceBetweenElectronTracks( + const reco::Track& track1, const reco::Track& track2, const edm::ESHandle& magneticField) const { + TwoTrackMinimumDistance ttmd; + FreeTrajectoryState fts1(GlobalPoint(track1.vx(), track1.vy(), track1.vz()), + GlobalVector(track1.px(), track1.py(), track1.pz()), + track1.charge(), + magneticField.product()); + FreeTrajectoryState fts2(GlobalPoint(track2.vx(), track2.vy(), track2.vz()), + GlobalVector(track2.px(), track2.py(), track2.pz()), + track2.charge(), + magneticField.product()); + bool status = ttmd.calculate(fts1, fts2); + if (!status) + return std::tuple(-999.f, status, GlobalPoint(-999.f, -999.f, -999.f)); + return std::make_tuple(ttmd.distance(), status, ttmd.crossingPoint()); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(ElectronVertexTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/JetImpactParameters.cc b/PhysicsTools/NanoAOD/plugins/JetImpactParameters.cc new file mode 100644 index 0000000000000..b5e081cc85ec3 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/JetImpactParameters.cc @@ -0,0 +1,140 @@ +/* +//////////////////////////////////////////////// + +Jet Impact parameter information for displaced tau collection : Pritam Palit, created on 01/09/2025 + +////////////////////////////////////////////////// + */ + +#include +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +void vector_test(std::vector& values) { + for (auto& value : values) { + if (std::isnan(value)) { + throw std::runtime_error("Jet IP output: NaN detected."); + } else if (std::isinf(value)) { + throw std::runtime_error("Jet IP output: Infinity detected."); + } else if (!std::isfinite(value)) { + throw std::runtime_error("Jet IP output: Non-standard value detected."); + } + } +} + +class JetImpactParameters : public edm::stream::EDProducer<> { +public: + explicit JetImpactParameters(const edm::ParameterSet&); + ~JetImpactParameters() override = default; + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + edm::EDGetTokenT jetsToken_; + edm::EDGetTokenT pfCandidatesToken_; + const double deltaRMax_; +}; + +JetImpactParameters::JetImpactParameters(const edm::ParameterSet& config) + : jetsToken_(consumes(config.getParameter("jets"))), + pfCandidatesToken_(consumes(config.getParameter("pfCandidates"))), + deltaRMax_(config.getParameter("deltaRMax")) { + produces>("jetDxy"); + produces>("jetDz"); + produces>("jetDxyError"); + produces>("jetDzError"); + produces>("jetCharge"); +} + +void JetImpactParameters::produce(edm::Event& event, const edm::EventSetup& setup) { + // Get jets and PFCandidates + auto jets = event.getHandle(jetsToken_); + + std::vector v_jetDxy(jets->size(), -1.0); + std::vector v_jetDz(jets->size(), -1.0); + std::vector v_jetDxyError(jets->size(), -1.0); + std::vector v_jetDzError(jets->size(), -1.0); + std::vector v_jetCharge(jets->size(), -1.0); + + // Loop over jets + for (size_t jetIndex = 0; jetIndex < jets->size(); ++jetIndex) { + const auto& jet = jets->at(jetIndex); + const auto& jetP4 = jet.polarP4(); + + // Find the leading charged PFCandidate within deltaR < 0.4 + const pat::PackedCandidate* leadingChargedPFCandidate = nullptr; + Float_t leadingPt = -1.0; + + // Loop over jet daughters + const size_t nDaughters = jet.numberOfDaughters(); + for (size_t i = 0; i < nDaughters; ++i) { + const auto& daughterPtr = jet.daughterPtr(i); + const auto* daughter = dynamic_cast(daughterPtr.get()); + + // Skip if not a charged candidate or doesn't have track details + if (!daughter || daughter->charge() == 0 || !daughter->hasTrackDetails()) + continue; + + Float_t deltaR = reco::deltaR(daughter->polarP4(), jetP4); + if (deltaR > deltaRMax_) + continue; + + if (daughter->pt() > leadingPt) { + leadingPt = daughter->pt(); + leadingChargedPFCandidate = daughter; + } + } + + if (leadingChargedPFCandidate) { + v_jetDxy.at(jetIndex) = leadingChargedPFCandidate->dxy(); + v_jetDz.at(jetIndex) = leadingChargedPFCandidate->dz(); + v_jetDxyError.at(jetIndex) = leadingChargedPFCandidate->dxyError(); + v_jetDzError.at(jetIndex) = leadingChargedPFCandidate->dzError(); + v_jetCharge.at(jetIndex) = leadingChargedPFCandidate->charge(); + } + } + + vector_test(v_jetDxy); + vector_test(v_jetDz); + vector_test(v_jetDxyError); + vector_test(v_jetDzError); + vector_test(v_jetCharge); + + auto vm_jetDxy = std::make_unique>(); + edm::ValueMap::Filler filler_jetDxy(*vm_jetDxy); + filler_jetDxy.insert(jets, v_jetDxy.begin(), v_jetDxy.end()); + filler_jetDxy.fill(); + event.put(std::move(vm_jetDxy), "jetDxy"); + + auto vm_jetDz = std::make_unique>(); + edm::ValueMap::Filler filler_jetDz(*vm_jetDz); + filler_jetDz.insert(jets, v_jetDz.begin(), v_jetDz.end()); + filler_jetDz.fill(); + event.put(std::move(vm_jetDz), "jetDz"); + + auto vm_jetDxyError = std::make_unique>(); + edm::ValueMap::Filler filler_jetDxyError(*vm_jetDxyError); + filler_jetDxyError.insert(jets, v_jetDxyError.begin(), v_jetDxyError.end()); + filler_jetDxyError.fill(); + event.put(std::move(vm_jetDxyError), "jetDxyError"); + + auto vm_jetDzError = std::make_unique>(); + edm::ValueMap::Filler filler_jetDzError(*vm_jetDzError); + filler_jetDzError.insert(jets, v_jetDzError.begin(), v_jetDzError.end()); + filler_jetDzError.fill(); + event.put(std::move(vm_jetDzError), "jetDzError"); + + auto vm_jetCharge = std::make_unique>(); + edm::ValueMap::Filler filler_jetCharge(*vm_jetCharge); + filler_jetCharge.insert(jets, v_jetCharge.begin(), v_jetCharge.end()); + filler_jetCharge.fill(); + event.put(std::move(vm_jetCharge), "jetCharge"); +} + +DEFINE_FWK_MODULE(JetImpactParameters); diff --git a/PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc b/PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc new file mode 100644 index 0000000000000..18d12db42cb74 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc @@ -0,0 +1,338 @@ +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonSelectors.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/NanoAOD/interface/FlatTable.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" + +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/IPTools/interface/IPTools.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +class MuonExtendedTableProducer : public edm::global::EDProducer<> { +public: + explicit MuonExtendedTableProducer(const edm::ParameterSet& iConfig) + : name_(iConfig.getParameter("name")), + rhoTag_(consumes(iConfig.getParameter("rho"))), + muonTag_(consumes>(iConfig.getParameter("muons"))), + vtxTag_(consumes(iConfig.getParameter("primaryVertex"))), + bsTag_(consumes(iConfig.getParameter("beamspot"))), + jetTag_(consumes>(iConfig.getParameter("jets"))), + jetFatTag_(consumes>(iConfig.getParameter("jetsFat"))), + jetSubTag_(consumes>(iConfig.getParameter("jetsSub"))), + transientTrackBuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) { + produces(); + } + + ~MuonExtendedTableProducer() override {}; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("rho")->setComment("input rho parameter"); + desc.add("muons")->setComment("input muon collection"); + desc.add("primaryVertex")->setComment("input primary vertex collection"); + desc.add("beamspot")->setComment("input beamspot collection"); + desc.add("jets")->setComment("input jet collection"); + desc.add("jetsFat")->setComment("input fat jet collection"); + desc.add("jetsSub")->setComment("input sub jet collection"); + desc.add("name")->setComment("name of the muon nanoaod::FlatTable we are extending"); + descriptions.add("muonTable", desc); + } + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + float getPFIso(const pat::Muon& muon) const; + int findMatchedJet(const reco::Candidate& lepton, const std::vector& jets) const; + void fillLeptonJetVariables(const reco::Muon* mu, + const std::vector& jets, + const reco::Vertex& vertex, + const double rho, + std::vector* jetIdx, + std::vector relIso0p4, + std::vector* jetPtRatio, + std::vector* jetPtRel, + std::vector* jetSelectedChargedMultiplicity) const; + + std::string name_; + edm::EDGetTokenT rhoTag_; + edm::EDGetTokenT> muonTag_; + edm::EDGetTokenT vtxTag_; + edm::EDGetTokenT bsTag_; + edm::EDGetTokenT> generalTrackTag_; + edm::EDGetTokenT> jetTag_; + edm::EDGetTokenT> jetFatTag_; + edm::EDGetTokenT> jetSubTag_; + edm::ESGetToken transientTrackBuilderToken_; +}; + +void MuonExtendedTableProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + const double& rho = iEvent.get(rhoTag_); + edm::Handle> muonHandle; + iEvent.getByToken(muonTag_, muonHandle); + const auto& muons = *muonHandle; + const reco::VertexCollection& primaryVertices = iEvent.get(vtxTag_); + const reco::BeamSpot& beamspots = iEvent.get(bsTag_); + const std::vector& jets = iEvent.get(jetTag_); + const std::vector& jetFat = iEvent.get(jetFatTag_); + const std::vector& jetSub = iEvent.get(jetSubTag_); + + const auto& pv = primaryVertices.at(0); + GlobalPoint primaryVertex(pv.x(), pv.y(), pv.z()); + + const auto& bs = beamspots.position(); + GlobalPoint beamSpot(bs.x(), bs.y(), bs.z()); + reco::Vertex beamSpotVertex(beamspots.position(), beamspots.covariance3D()); + + const TransientTrackBuilder& builder = iSetup.getData(transientTrackBuilderToken_); + + unsigned int nMuons = muons.size(); + + std::vector idx, charge, trkPt, trkPtErr; + std::vector numberInnerHitsMissing; + std::vector relIso0p4; + std::vector jetPtRatio, jetPtRel; + std::vector jetIdx; + std::vector jetFatIdx, jetSubIdx; + std::vector jetSelectedChargedMultiplicity; + + std::vector dzPV, dzPVErr, dxyPVTraj, dxyPVTrajErr, dxyPVSigned, dxyPVSignedErr, ip3DPVSigned, ip3DPVSignedErr; + + std::vector trkNumPlanes, trkNumHits, trkNumDTHits, trkNumCSCHits, trkNumPixelHits(nMuons, -1), + trkNumTrkLayers(nMuons, -1), normChi2; + std::vector outerEta(nMuons, -5), outerPhi(nMuons, -5); + + for (unsigned int i = 0; i < nMuons; i++) { + const pat::Muon& muon = muons[i]; + const pat::MuonRef muonRef(muonHandle, i); + + reco::TrackRef trackRef; + + if (muon.isGlobalMuon()) { + trackRef = muon.combinedMuon(); + } else if (muon.isStandAloneMuon()) { + trackRef = muon.standAloneMuon(); + } else { + trackRef = muon.tunePMuonBestTrack(); + } + + idx.push_back(i); + + const auto& track = trackRef.get(); + reco::TransientTrack transientTrack = builder.build(track); + + charge.push_back(muon.charge()); + + trkPt.push_back(track->pt()); + trkPtErr.push_back(track->ptError()); + + relIso0p4.push_back(getPFIso(muon)); + + fillLeptonJetVariables( + &muon, jets, pv, rho, &jetIdx, relIso0p4, &jetPtRatio, &jetPtRel, &jetSelectedChargedMultiplicity); + + const reco::Candidate* mu_cand = dynamic_cast(&muon); + jetFatIdx.push_back(findMatchedJet(*mu_cand, jetFat)); + jetSubIdx.push_back(findMatchedJet(*mu_cand, jetSub)); + + numberInnerHitsMissing.push_back(!muon.innerTrack().isNull() ? muon.innerTrack()->hitPattern().numberOfLostHits( + reco::HitPattern::MISSING_INNER_HITS) + : 0); + + dzPV.push_back(track->dz(pv.position())); + dzPVErr.push_back(std::hypot(track->dzError(), pv.zError())); + TrajectoryStateClosestToPoint trajectoryPV = transientTrack.trajectoryStateClosestToPoint(primaryVertex); + dxyPVTraj.push_back(trajectoryPV.perigeeParameters().transverseImpactParameter()); + dxyPVTrajErr.push_back(trajectoryPV.perigeeError().transverseImpactParameterError()); + GlobalVector muonRefTrackDir(muon.px(), muon.py(), muon.pz()); + dxyPVSigned.push_back(IPTools::signedTransverseImpactParameter(transientTrack, muonRefTrackDir, pv).second.value()); + dxyPVSignedErr.push_back( + IPTools::signedTransverseImpactParameter(transientTrack, muonRefTrackDir, pv).second.error()); + + ip3DPVSigned.push_back( + IPTools::signedImpactParameter3D(transientTrack, muonRefTrackDir, beamSpotVertex).second.value()); + ip3DPVSignedErr.push_back( + IPTools::signedImpactParameter3D(transientTrack, muonRefTrackDir, beamSpotVertex).second.error()); + + trkNumPlanes.push_back(track->hitPattern().muonStationsWithValidHits()); + trkNumHits.push_back(track->hitPattern().numberOfValidMuonHits()); + trkNumDTHits.push_back(track->hitPattern().numberOfValidMuonDTHits()); + trkNumCSCHits.push_back(track->hitPattern().numberOfValidMuonCSCHits()); + + normChi2.push_back(track->normalizedChi2()); + + // Cannot get outer track for tracker muons + if (track->extra().isNonnull() && track->extra().isAvailable() && track->outerOk()) { + outerEta[i] = track->outerEta(); + outerPhi[i] = track->outerPhi(); + } + } + + auto tab = std::make_unique(nMuons, name_, false, true); + tab->addColumn("idx", idx, "EXOnanoAOD muon index"); + + tab->addColumn("trkPt", trkPt, ""); + tab->addColumn("trkPtErr", trkPtErr, ""); + + tab->addColumn("relIso0p4", relIso0p4, ""); + tab->addColumn("jetPtRatio", jetPtRatio, ""); + tab->addColumn("jetPtRel", jetPtRel, ""); + tab->addColumn("jetSelectedChargedMultiplicity", jetSelectedChargedMultiplicity, ""); + tab->addColumn("jetIdx", jetIdx, ""); + tab->addColumn("jetFatIdx", jetFatIdx, ""); + tab->addColumn("jetSubIdx", jetSubIdx, ""); + + tab->addColumn("numberInnerHitsMissing", numberInnerHitsMissing, ""); + + tab->addColumn("dzPV", dzPV, ""); + tab->addColumn("dzPVErr", dzPVErr, ""); + tab->addColumn("dxyPVTraj", dxyPVTraj, ""); + tab->addColumn("dxyPVTrajErr", dxyPVTrajErr, ""); + tab->addColumn("dxyPVSigned", dxyPVSigned, ""); + tab->addColumn("dxyPVSignedErr", dxyPVSignedErr, ""); + tab->addColumn("ip3DPVSigned", ip3DPVSigned, ""); + tab->addColumn("ip3DPVSignedErr", ip3DPVSignedErr, ""); + + tab->addColumn("trkNumPlanes", trkNumPlanes, ""); + tab->addColumn("trkNumHits", trkNumHits, ""); + tab->addColumn("trkNumDTHits", trkNumDTHits, ""); + tab->addColumn("trkNumCSCHits", trkNumCSCHits, ""); + tab->addColumn("normChi2", normChi2, ""); + tab->addColumn("trkNumPixelHits", trkNumPixelHits, ""); + tab->addColumn("trkNumTrkLayers", trkNumTrkLayers, ""); + + tab->addColumn("outerEta", outerEta, ""); + tab->addColumn("outerPhi", outerPhi, ""); + + iEvent.put(std::move(tab)); +} + +float MuonExtendedTableProducer::getPFIso(const pat::Muon& muon) const { + return (muon.pfIsolationR04().sumChargedHadronPt + + std::max(0., + muon.pfIsolationR04().sumNeutralHadronEt + muon.pfIsolationR04().sumPhotonEt - + 0.5 * muon.pfIsolationR04().sumPUPt)) / + muon.pt(); +} + +template +bool isSourceCandidatePtrMatch(const T1& lhs, const T2& rhs) { + for (size_t lhsIndex = 0; lhsIndex < lhs.numberOfSourceCandidatePtrs(); ++lhsIndex) { + auto lhsSourcePtr = lhs.sourceCandidatePtr(lhsIndex); + for (size_t rhsIndex = 0; rhsIndex < rhs.numberOfSourceCandidatePtrs(); ++rhsIndex) { + auto rhsSourcePtr = rhs.sourceCandidatePtr(rhsIndex); + if (lhsSourcePtr == rhsSourcePtr) { + return true; + } + } + } + + return false; +} + +int MuonExtendedTableProducer::findMatchedJet(const reco::Candidate& lepton, const std::vector& jets) const { + int iJet = -1; + + unsigned int nJets = jets.size(); + + for (unsigned int i = 0; i < nJets; i++) { + const pat::Jet& jet = jets[i]; + if (isSourceCandidatePtrMatch(lepton, jet)) { + return i; + } + } + + return iJet; +} + +void MuonExtendedTableProducer::fillLeptonJetVariables(const reco::Muon* mu, + const std::vector& jets, + const reco::Vertex& vertex, + const double rho, + std::vector* jetIdx, + std::vector relIso0p4, + std::vector* jetPtRatio, + std::vector* jetPtRel, + std::vector* jetSelectedChargedMultiplicity) const { + const reco::Candidate* cand = dynamic_cast(mu); + int matchedJetIdx = findMatchedJet(*cand, jets); + + jetIdx->push_back(matchedJetIdx); + + if (matchedJetIdx < 0) { + float ptRatio = (1. / (1. + relIso0p4.back())); + jetPtRatio->push_back(ptRatio); + jetPtRel->push_back(0); + jetSelectedChargedMultiplicity->push_back(0); + } else { + const pat::Jet& jet = jets[matchedJetIdx]; + auto rawJetP4 = jet.correctedP4("Uncorrected"); + auto leptonP4 = cand->p4(); + + bool leptonEqualsJet = ((rawJetP4 - leptonP4).P() < 1e-4); + + if (leptonEqualsJet) { + jetPtRatio->push_back(1); + jetPtRel->push_back(0); + jetSelectedChargedMultiplicity->push_back(0); + } else { + auto L1JetP4 = jet.correctedP4("L1FastJet"); + double L2L3JEC = jet.pt() / L1JetP4.pt(); + auto lepAwareJetP4 = (L1JetP4 - leptonP4) * L2L3JEC + leptonP4; + + float ptRatio = cand->pt() / lepAwareJetP4.pt(); + float ptRel = leptonP4.Vect().Cross((lepAwareJetP4 - leptonP4).Vect().Unit()).R(); + jetPtRatio->push_back(ptRatio); + jetPtRel->push_back(ptRel); + jetSelectedChargedMultiplicity->push_back(0); + + for (const auto& daughterPtr : jet.daughterPtrVector()) { + const pat::PackedCandidate& daughter = *((const pat::PackedCandidate*)daughterPtr.get()); + + if (daughter.charge() == 0) + continue; + if (daughter.fromPV() < 2) + continue; + if (reco::deltaR(daughter, *cand) > 0.4) + continue; + if (!daughter.hasTrackDetails()) + continue; + + auto daughterTrack = daughter.pseudoTrack(); + + if (daughterTrack.pt() <= 1) + continue; + if (daughterTrack.hitPattern().numberOfValidHits() < 8) + continue; + if (daughterTrack.hitPattern().numberOfValidPixelHits() < 2) + continue; + if (daughterTrack.normalizedChi2() >= 5) + continue; + if (std::abs(daughterTrack.dz(vertex.position())) >= 17) + continue; + if (std::abs(daughterTrack.dxy(vertex.position())) >= 0.2) + continue; + ++jetSelectedChargedMultiplicity->back(); + } + } + } +} + +#include "FWCore/Framework/interface/MakerMacros.h" +//define this as a plug-in +DEFINE_FWK_MODULE(MuonExtendedTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/MuonVertexTableProducer.cc b/PhysicsTools/NanoAOD/plugins/MuonVertexTableProducer.cc new file mode 100644 index 0000000000000..c945e9e733b7f --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/MuonVertexTableProducer.cc @@ -0,0 +1,874 @@ +#include "FWCore/Framework/interface/global/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/transform.h" +#include "DataFormats/NanoAOD/interface/FlatTable.h" + +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "TLorentzVector.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/PatCandidates/interface/Muon.h" + +#include "DataFormats/Math/interface/deltaR.h" + +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h" + +#include "TrackingTools/Records/interface/TrackingComponentsRecord.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" +#include "TrackingTools/GeomPropagators/interface/Propagator.h" + +#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h" +#include "PhysicsTools/RecoUtils/interface/CheckHitPattern.h" +#include "DataFormats/TrackerCommon/interface/TrackerTopology.h" +#include "Geometry/Records/interface/TrackerTopologyRcd.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" + +#include "TrackingTools/IPTools/interface/IPTools.h" + +#include +#include + +class MuonVertexTableProducer : public edm::global::EDProducer<> { +public: + explicit MuonVertexTableProducer(const edm::ParameterSet& iConfig) + : dsaMuonTag_(consumes>(iConfig.getParameter("dsaMuons"))), + patMuonTag_(consumes>(iConfig.getParameter("patMuons"))), + bsTag_(consumes(iConfig.getParameter("beamspot"))), + pvTag_(consumes(iConfig.getParameter("primaryVertex"))), + generalTrackTag_(consumes>(iConfig.getParameter("generalTracks"))), + propagatorToken_(esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))), + magneticFieldToken_(esConsumes()), + tkerGeomToken_(esConsumes()), + tkerTopoToken_(esConsumes()), + transientTrackBuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) { + produces("PatMuonVertex"); + produces("PatDSAMuonVertex"); + produces("DSAMuonVertex"); + produces("PatMuonVertexRefittedTracks"); + produces("PatDSAMuonVertexRefittedTracks"); + produces("DSAMuonVertexRefittedTracks"); + } + + ~MuonVertexTableProducer() override {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("patMuons")->setComment("input pat muon collection"); + desc.add("dsaMuons")->setComment("input displaced standalone muon collection"); + desc.add("beamspot")->setComment("input beamspot collection"); + desc.add("primaryVertex")->setComment("input primaryVertex collection"); + desc.add("generalTracks")->setComment("input generalTracks collection"); + descriptions.add("muonVertexTables", desc); + } + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + std::pair getVxy(const reco::Vertex muonVertex) const; + std::pair getVxyz(const reco::Vertex muonVertex) const; + + template + float getDisplacedTrackerIsolation(const std::vector& generalTracks, + const MuonType1& muon_1, + const reco::Vertex muonVertex, + const reco::BeamSpot& beamspot, + const MuonType2* muon_2 = nullptr, + float maxDR = 0.3, + float minDR = 0.01, + float maxDz = 0.5, + float maxDxy = 0.2) const; + + template + float getProximityDeltaR(const MuonType& track, + const MuonType& trackRef, + const edm::ESHandle& magneticField, + const edm::ESHandle& propagator) const; + + std::tuple getDistanceBetweenMuonTracks( + const reco::Track& track1, const reco::Track& track2, const edm::ESHandle& magneticField) const; + + const edm::EDGetTokenT> dsaMuonTag_; + const edm::EDGetTokenT> patMuonTag_; + const edm::EDGetTokenT bsTag_; + const edm::EDGetTokenT pvTag_; + const edm::EDGetTokenT> generalTrackTag_; + const edm::ESGetToken propagatorToken_; + const edm::ESGetToken magneticFieldToken_; + const edm::ESGetToken tkerGeomToken_; + const edm::ESGetToken tkerTopoToken_; + const edm::ESGetToken transientTrackBuilderToken_; +}; + +void MuonVertexTableProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + const std::vector& dsaMuons = iEvent.get(dsaMuonTag_); + const std::vector& patMuons = iEvent.get(patMuonTag_); + + const reco::BeamSpot& beamSpotInput = iEvent.get(bsTag_); + const auto& bs = beamSpotInput.position(); + GlobalPoint beamSpot(bs.x(), bs.y(), bs.z()); + reco::Vertex beamSpotVertex(beamSpotInput.position(), beamSpotInput.covariance3D()); + + const reco::VertexCollection& primaryVertices = iEvent.get(pvTag_); + const auto& pv = primaryVertices.at(0); + GlobalPoint primaryVertex(pv.x(), pv.y(), pv.z()); + + const std::vector& generalTracks = iEvent.get(generalTrackTag_); + + auto const& propagator = &iSetup.getData(propagatorToken_); + auto const& magneticField = &iSetup.getData(magneticFieldToken_); + + auto const& tkerGeom = &iSetup.getData(tkerGeomToken_); + auto const& tkerTopo = &iSetup.getData(tkerTopoToken_); + const TransientTrackBuilder& builder = iSetup.getData(transientTrackBuilderToken_); + + KalmanVertexFitter vertexFitter(true, true); + + int nPatPatVertices = 0; + int nPatDSAVertices = 0; + int nDSADSAVertices = 0; + int ppRefittedTrackIdx_counter = 0; + int pdRefittedTrackIdx_counter = 0; + int ddRefittedTrackIdx_counter = 0; + + std::map> vertexIsValid; + std::map> vxy, vxyz, vx, vy, vz, t, vxySigma, vxyzSigma, vxErr, vyErr, vzErr, tErr; + std::map> chi2, ndof, normChi2, dR, originalMuonIdx1, originalMuonIdx2, + refittedTrackIdx1, refittedTrackIdx2, isDSAMuon1, isDSAMuon2; + std::map> displacedTrackIso03Dimuon1, displacedTrackIso03Dimuon2, + displacedTrackIso04Dimuon1, displacedTrackIso04Dimuon2; + std::map> displacedTrackIso03Muon1, displacedTrackIso03Muon2, + displacedTrackIso04Muon1, displacedTrackIso04Muon2, proxDeltaR; + std::map> DCA, DCAstatus, DCAx, DCAy, DCAz; + std::map> hitsInFrontOfVert1, missHitsAfterVert1, hitsInFrontOfVert2, + missHitsAfterVert2; + + std::map> refittedTrackIdx, refittedTrackIsDSAMuon, refittedTrackOriginalIdx, + refittedTrackPt, refittedTrackPtErr, refittedTrackPx, refittedTrackPy, refittedTrackPz; + std::map> refittedTrackEta, refittedTrackEtaErr, refittedTrackPhi, + refittedTrackPhiErr, refittedTrackCharge, refittedTrackNormChi2, refittedTrackNdof, refittedTrackChi2; + std::map> refittedTrackDzPV, refittedTrackDzPVErr, refittedTrackDxyPVTraj, + refittedTrackDxyPVTrajErr, refittedTrackDxyPVSigned, refittedTrackDxyPVSignedErr, refittedTrackIp3DPVSigned, + refittedTrackIp3DPVSignedErr; + std::map> refittedTrackIso03Dimuon1, refittedTrackIso03Dimuon2, + refittedTrackIso04Dimuon1, refittedTrackIso04Dimuon2, refittedTrackIso03Muon1, refittedTrackIso03Muon2, + refittedTrackIso04Muon1, refittedTrackIso04Muon2; + + // pat muons + for (size_t i = 0; i < patMuons.size(); i++) { + const pat::Muon& muon_i = patMuons.at(i); + + reco::TrackRef trackRef_i; + if (muon_i.isGlobalMuon()) + trackRef_i = muon_i.combinedMuon(); + else if (muon_i.isStandAloneMuon()) + trackRef_i = muon_i.standAloneMuon(); + else + trackRef_i = muon_i.tunePMuonBestTrack(); + + const auto& muonTrack_i = trackRef_i.get(); + reco::TransientTrack muonTransientTrack_i = builder.build(muonTrack_i); + + // pat-pat muon vertex + for (size_t j = i + 1; j < patMuons.size(); j++) { + const pat::Muon& muon_j = patMuons.at(j); + + reco::TrackRef trackRef_j; + if (muon_j.isGlobalMuon()) + trackRef_j = muon_j.combinedMuon(); + else if (muon_j.isStandAloneMuon()) + trackRef_j = muon_j.standAloneMuon(); + else + trackRef_j = muon_j.tunePMuonBestTrack(); + + const auto& muonTrack_j = trackRef_j.get(); + reco::TransientTrack muonTransientTrack_j = builder.build(muonTrack_j); + + std::vector muonTransientTracks{}; + muonTransientTracks.push_back(muonTransientTrack_i); + muonTransientTracks.push_back(muonTransientTrack_j); + + TransientVertex transientMuonVertex = vertexFitter.vertex(muonTransientTracks); + reco::Vertex muonVertex = reco::Vertex(transientMuonVertex); + + if (!transientMuonVertex.isValid()) + continue; + std::tuple distanceTuple = + getDistanceBetweenMuonTracks(*muonTrack_i, *muonTrack_j, magneticField); + // if dca status is good but dca is more than 15 cm + if (std::get<1>(distanceTuple) && std::get<0>(distanceTuple) > 15) + continue; + + nPatPatVertices++; + + vertexIsValid["PATPAT"].push_back(transientMuonVertex.isValid()); + std::pair vxy_ = getVxy(muonVertex); + vxy["PATPAT"].push_back(vxy_.first); + vxySigma["PATPAT"].push_back(vxy_.second); + std::pair vxyz_ = getVxyz(muonVertex); + vxyz["PATPAT"].push_back(vxyz_.first); + vxyzSigma["PATPAT"].push_back(vxyz_.second); + chi2["PATPAT"].push_back(muonVertex.chi2()); + ndof["PATPAT"].push_back(muonVertex.ndof()); + normChi2["PATPAT"].push_back(muonVertex.normalizedChi2()); + vx["PATPAT"].push_back(muonVertex.x()); + vy["PATPAT"].push_back(muonVertex.y()); + vz["PATPAT"].push_back(muonVertex.z()); + t["PATPAT"].push_back(muonVertex.t()); + vxErr["PATPAT"].push_back(muonVertex.xError()); + vyErr["PATPAT"].push_back(muonVertex.yError()); + vzErr["PATPAT"].push_back(muonVertex.zError()); + tErr["PATPAT"].push_back(muonVertex.tError()); + + dR["PATPAT"].push_back(reco::deltaR(muon_i, muon_j)); + originalMuonIdx1["PATPAT"].push_back(i); + originalMuonIdx2["PATPAT"].push_back(j); + isDSAMuon1["PATPAT"].push_back(0); + isDSAMuon2["PATPAT"].push_back(0); + + displacedTrackIso03Dimuon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, &muon_j, 0.3)); + displacedTrackIso04Dimuon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, &muon_j, 0.4)); + displacedTrackIso03Dimuon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_j, muonVertex, beamSpotInput, &muon_i, 0.3)); + displacedTrackIso04Dimuon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_j, muonVertex, beamSpotInput, &muon_i, 0.4)); + displacedTrackIso03Muon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, nullptr, 0.3)); + displacedTrackIso04Muon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, nullptr, 0.4)); + displacedTrackIso03Muon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_j, muonVertex, beamSpotInput, nullptr, 0.3)); + displacedTrackIso04Muon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_j, muonVertex, beamSpotInput, nullptr, 0.4)); + + // PAT muons cannot be only tracker muons for proximity deltaR calculations + if ((muon_i.isGlobalMuon() || muon_i.isStandAloneMuon()) && + (muon_j.isGlobalMuon() || muon_j.isStandAloneMuon())) { + proxDeltaR["PATPAT"].push_back(getProximityDeltaR(*muonTrack_i, *muonTrack_j, magneticField, propagator)); + } else + proxDeltaR["PATPAT"].push_back(-1); + + DCA["PATPAT"].push_back(std::get<0>(distanceTuple)); + DCAstatus["PATPAT"].push_back(std::get<1>(distanceTuple)); + DCAx["PATPAT"].push_back(std::get<2>(distanceTuple).x()); + DCAy["PATPAT"].push_back(std::get<2>(distanceTuple).y()); + DCAz["PATPAT"].push_back(std::get<2>(distanceTuple).z()); + + CheckHitPattern checkHitPattern; + checkHitPattern.init(tkerTopo, *tkerGeom, builder); + if (muon_i.isTrackerMuon()) { + CheckHitPattern::Result hitPattern_i = checkHitPattern(*muonTrack_i, transientMuonVertex.vertexState()); + hitsInFrontOfVert1["PATPAT"].push_back(hitPattern_i.hitsInFrontOfVert); + missHitsAfterVert1["PATPAT"].push_back(hitPattern_i.missHitsAfterVert); + } else { + hitsInFrontOfVert1["PATPAT"].push_back(-1); + missHitsAfterVert1["PATPAT"].push_back(-1); + } + if (muon_j.isTrackerMuon()) { + CheckHitPattern::Result hitPattern_j = checkHitPattern(*muonTrack_j, transientMuonVertex.vertexState()); + hitsInFrontOfVert2["PATPAT"].push_back(hitPattern_j.hitsInFrontOfVert); + missHitsAfterVert2["PATPAT"].push_back(hitPattern_j.missHitsAfterVert); + } else { + hitsInFrontOfVert2["PATPAT"].push_back(-1); + missHitsAfterVert2["PATPAT"].push_back(-1); + } + + reco::TransientTrack refittedTrack_i = transientMuonVertex.refittedTrack(muonTransientTrack_i); + reco::TransientTrack refittedTrack_j = transientMuonVertex.refittedTrack(muonTransientTrack_j); + std::vector refittedTracks = {refittedTrack_i, refittedTrack_j}; + for (size_t k = 0; k < refittedTracks.size(); ++k) { + const auto& refittedTrack = refittedTracks[k]; + refittedTrackIsDSAMuon["PATPAT"].push_back(0); + refittedTrackOriginalIdx["PATPAT"].push_back(i); + refittedTrackPt["PATPAT"].push_back(refittedTrack.track().pt()); + refittedTrackPtErr["PATPAT"].push_back(refittedTrack.track().ptError()); + refittedTrackPx["PATPAT"].push_back(refittedTrack.track().px()); + refittedTrackPy["PATPAT"].push_back(refittedTrack.track().py()); + refittedTrackPz["PATPAT"].push_back(refittedTrack.track().pz()); + refittedTrackEta["PATPAT"].push_back(refittedTrack.track().eta()); + refittedTrackEtaErr["PATPAT"].push_back(refittedTrack.track().etaError()); + refittedTrackPhi["PATPAT"].push_back(refittedTrack.track().phi()); + refittedTrackPhiErr["PATPAT"].push_back(refittedTrack.track().phiError()); + refittedTrackCharge["PATPAT"].push_back(refittedTrack.track().charge()); + refittedTrackNormChi2["PATPAT"].push_back(refittedTrack.normalizedChi2()); + refittedTrackNdof["PATPAT"].push_back(refittedTrack.ndof()); + refittedTrackChi2["PATPAT"].push_back(refittedTrack.chi2()); + + refittedTrackDzPV["PATPAT"].push_back(refittedTrack.track().dz(pv.position())); + refittedTrackDzPVErr["PATPAT"].push_back(std::hypot(refittedTrack.track().dzError(), pv.zError())); + TrajectoryStateClosestToPoint trajectoryPV_i = refittedTrack.trajectoryStateClosestToPoint(primaryVertex); + refittedTrackDxyPVTraj["PATPAT"].push_back(trajectoryPV_i.perigeeParameters().transverseImpactParameter()); + refittedTrackDxyPVTrajErr["PATPAT"].push_back(trajectoryPV_i.perigeeError().transverseImpactParameterError()); + GlobalVector muonRefTrackDir_i( + refittedTrack.track().px(), refittedTrack.track().py(), refittedTrack.track().pz()); + refittedTrackDxyPVSigned["PATPAT"].push_back( + IPTools::signedTransverseImpactParameter(refittedTrack, muonRefTrackDir_i, pv).second.value()); + refittedTrackDxyPVSignedErr["PATPAT"].push_back( + IPTools::signedTransverseImpactParameter(refittedTrack, muonRefTrackDir_i, pv).second.error()); + refittedTrackIp3DPVSigned["PATPAT"].push_back( + IPTools::signedImpactParameter3D(refittedTrack, muonRefTrackDir_i, pv).second.value()); + refittedTrackIp3DPVSignedErr["PATPAT"].push_back( + IPTools::signedImpactParameter3D(refittedTrack, muonRefTrackDir_i, pv).second.error()); + + refittedTrackIdx["PATPAT"].push_back(ppRefittedTrackIdx_counter); + if (k == 0) + refittedTrackIdx1["PATPAT"].push_back(ppRefittedTrackIdx_counter); + if (k == 1) + refittedTrackIdx2["PATPAT"].push_back(ppRefittedTrackIdx_counter); + ppRefittedTrackIdx_counter++; + } + + refittedTrackIso03Dimuon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, &refittedTrack_j.track(), 0.3)); + refittedTrackIso04Dimuon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, &refittedTrack_j.track(), 0.4)); + refittedTrackIso03Dimuon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, &refittedTrack_i.track(), 0.3)); + refittedTrackIso04Dimuon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, &refittedTrack_i.track(), 0.4)); + refittedTrackIso03Muon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, nullptr, 0.3)); + refittedTrackIso04Muon1["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, nullptr, 0.4)); + refittedTrackIso03Muon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, nullptr, 0.3)); + refittedTrackIso04Muon2["PATPAT"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, nullptr, 0.4)); + } + + // pat-dsa muon vertex + for (size_t j = 0; j < dsaMuons.size(); j++) { + const auto& muonTrack_j = dsaMuons.at(j); + reco::TransientTrack muonTransientTrack_j = builder.build(muonTrack_j); + + std::vector muonTransientTracks{}; + muonTransientTracks.push_back(muonTransientTrack_i); + muonTransientTracks.push_back(muonTransientTrack_j); + + TransientVertex transientMuonVertex = vertexFitter.vertex(muonTransientTracks); + reco::Vertex muonVertex = reco::Vertex(transientMuonVertex); + + if (!transientMuonVertex.isValid()) + continue; + std::tuple distanceTuple = + getDistanceBetweenMuonTracks(*muonTrack_i, muonTrack_j, magneticField); + // if dca status is good but dca is more than 15 cm + if (std::get<1>(distanceTuple) && std::get<0>(distanceTuple) > 15) + continue; + + nPatDSAVertices++; + + vertexIsValid["PATDSA"].push_back(transientMuonVertex.isValid()); + std::pair vxy_ = getVxy(muonVertex); + vxy["PATDSA"].push_back(vxy_.first); + vxySigma["PATDSA"].push_back(vxy_.second); + std::pair vxyz_ = getVxyz(muonVertex); + vxyz["PATDSA"].push_back(vxyz_.first); + vxyzSigma["PATDSA"].push_back(vxyz_.second); + chi2["PATDSA"].push_back(muonVertex.chi2()); + ndof["PATDSA"].push_back(muonVertex.ndof()); + normChi2["PATDSA"].push_back(muonVertex.normalizedChi2()); + vx["PATDSA"].push_back(muonVertex.x()); + vy["PATDSA"].push_back(muonVertex.y()); + vz["PATDSA"].push_back(muonVertex.z()); + t["PATDSA"].push_back(muonVertex.t()); + vxErr["PATDSA"].push_back(muonVertex.xError()); + vyErr["PATDSA"].push_back(muonVertex.yError()); + vzErr["PATDSA"].push_back(muonVertex.zError()); + tErr["PATDSA"].push_back(muonVertex.tError()); + + dR["PATDSA"].push_back(reco::deltaR(muon_i.eta(), muonTrack_j.eta(), muon_i.phi(), muonTrack_j.phi())); + originalMuonIdx1["PATDSA"].push_back(i); + originalMuonIdx2["PATDSA"].push_back(j); + isDSAMuon1["PATDSA"].push_back(0); + isDSAMuon2["PATDSA"].push_back(1); + + displacedTrackIso03Dimuon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, &muonTrack_j, 0.3)); + displacedTrackIso04Dimuon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, &muonTrack_j, 0.4)); + displacedTrackIso03Dimuon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, &muon_i, 0.3)); + displacedTrackIso04Dimuon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, &muon_i, 0.4)); + displacedTrackIso03Muon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, nullptr, 0.3)); + displacedTrackIso04Muon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muon_i, muonVertex, beamSpotInput, nullptr, 0.4)); + displacedTrackIso03Muon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, nullptr, 0.3)); + displacedTrackIso04Muon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, nullptr, 0.4)); + + // PAT muons cannot be only tracker muons for proximity deltaR calculations + if (muon_i.isGlobalMuon() || muon_i.isStandAloneMuon()) { + proxDeltaR["PATDSA"].push_back(getProximityDeltaR(*muonTrack_i, muonTrack_j, magneticField, propagator)); + } else + proxDeltaR["PATDSA"].push_back(-1); + + DCA["PATDSA"].push_back(std::get<0>(distanceTuple)); + DCAstatus["PATDSA"].push_back(std::get<1>(distanceTuple)); + DCAx["PATDSA"].push_back(std::get<2>(distanceTuple).x()); + DCAy["PATDSA"].push_back(std::get<2>(distanceTuple).y()); + DCAz["PATDSA"].push_back(std::get<2>(distanceTuple).z()); + + CheckHitPattern checkHitPattern; + checkHitPattern.init(tkerTopo, *tkerGeom, builder); + if (muon_i.isTrackerMuon()) { + CheckHitPattern::Result hitPattern_i = checkHitPattern(*muonTrack_i, transientMuonVertex.vertexState()); + hitsInFrontOfVert1["PATDSA"].push_back(hitPattern_i.hitsInFrontOfVert); + missHitsAfterVert1["PATDSA"].push_back(hitPattern_i.missHitsAfterVert); + } else { + hitsInFrontOfVert1["PATDSA"].push_back(-1); + missHitsAfterVert1["PATDSA"].push_back(-1); + } + hitsInFrontOfVert2["PATDSA"].push_back(-1); + missHitsAfterVert2["PATDSA"].push_back(-1); + + reco::TransientTrack refittedTrack_i = transientMuonVertex.refittedTrack(muonTransientTrack_i); + reco::TransientTrack refittedTrack_j = transientMuonVertex.refittedTrack(muonTransientTrack_j); + std::vector refittedTracks = {refittedTrack_i, refittedTrack_j}; + for (size_t k = 0; k < refittedTracks.size(); ++k) { + const auto& refittedTrack = refittedTracks[k]; + refittedTrackIsDSAMuon["PATDSA"].push_back(0); + refittedTrackOriginalIdx["PATDSA"].push_back(i); + refittedTrackPt["PATDSA"].push_back(refittedTrack.track().pt()); + refittedTrackPtErr["PATDSA"].push_back(refittedTrack.track().ptError()); + refittedTrackPx["PATDSA"].push_back(refittedTrack.track().px()); + refittedTrackPy["PATDSA"].push_back(refittedTrack.track().py()); + refittedTrackPz["PATDSA"].push_back(refittedTrack.track().pz()); + refittedTrackEta["PATDSA"].push_back(refittedTrack.track().eta()); + refittedTrackEtaErr["PATDSA"].push_back(refittedTrack.track().etaError()); + refittedTrackPhi["PATDSA"].push_back(refittedTrack.track().phi()); + refittedTrackPhiErr["PATDSA"].push_back(refittedTrack.track().phiError()); + refittedTrackCharge["PATDSA"].push_back(refittedTrack.track().charge()); + refittedTrackNormChi2["PATDSA"].push_back(refittedTrack.normalizedChi2()); + refittedTrackNdof["PATDSA"].push_back(refittedTrack.ndof()); + refittedTrackChi2["PATDSA"].push_back(refittedTrack.chi2()); + + refittedTrackDzPV["PATDSA"].push_back(refittedTrack.track().dz(pv.position())); + refittedTrackDzPVErr["PATDSA"].push_back(std::hypot(refittedTrack.track().dzError(), pv.zError())); + TrajectoryStateClosestToPoint trajectoryPV_i = refittedTrack.trajectoryStateClosestToPoint(primaryVertex); + refittedTrackDxyPVTraj["PATDSA"].push_back(trajectoryPV_i.perigeeParameters().transverseImpactParameter()); + refittedTrackDxyPVTrajErr["PATDSA"].push_back(trajectoryPV_i.perigeeError().transverseImpactParameterError()); + GlobalVector muonRefTrackDir_i( + refittedTrack.track().px(), refittedTrack.track().py(), refittedTrack.track().pz()); + refittedTrackDxyPVSigned["PATDSA"].push_back( + IPTools::signedTransverseImpactParameter(refittedTrack, muonRefTrackDir_i, pv).second.value()); + refittedTrackDxyPVSignedErr["PATDSA"].push_back( + IPTools::signedTransverseImpactParameter(refittedTrack, muonRefTrackDir_i, pv).second.error()); + refittedTrackIp3DPVSigned["PATDSA"].push_back( + IPTools::signedImpactParameter3D(refittedTrack, muonRefTrackDir_i, pv).second.value()); + refittedTrackIp3DPVSignedErr["PATDSA"].push_back( + IPTools::signedImpactParameter3D(refittedTrack, muonRefTrackDir_i, pv).second.error()); + + refittedTrackIdx["PATDSA"].push_back(pdRefittedTrackIdx_counter); + if (k == 0) + refittedTrackIdx1["PATDSA"].push_back(pdRefittedTrackIdx_counter); + if (k == 1) + refittedTrackIdx2["PATDSA"].push_back(pdRefittedTrackIdx_counter); + pdRefittedTrackIdx_counter++; + } + + refittedTrackIso03Dimuon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, &refittedTrack_j.track(), 0.3)); + refittedTrackIso04Dimuon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, &refittedTrack_j.track(), 0.4)); + refittedTrackIso03Dimuon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, &refittedTrack_i.track(), 0.3)); + refittedTrackIso04Dimuon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, &refittedTrack_i.track(), 0.4)); + refittedTrackIso03Muon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, nullptr, 0.3)); + refittedTrackIso04Muon1["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, nullptr, 0.4)); + refittedTrackIso03Muon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, nullptr, 0.3)); + refittedTrackIso04Muon2["PATDSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, nullptr, 0.4)); + } + } + + // dsa muons + for (size_t i = 0; i < dsaMuons.size(); i++) { + const auto& muonTrack_i = dsaMuons.at(i); + reco::TransientTrack muonTransientTrack_i = builder.build(muonTrack_i); + + //dsa-dsa muon vertex + for (size_t j = i + 1; j < dsaMuons.size(); j++) { + const auto& muonTrack_j = dsaMuons.at(j); + reco::TransientTrack muonTransientTrack_j = builder.build(muonTrack_j); + + std::vector muonTransientTracks{}; + muonTransientTracks.push_back(muonTransientTrack_i); + muonTransientTracks.push_back(muonTransientTrack_j); + + TransientVertex transientMuonVertex = vertexFitter.vertex(muonTransientTracks); + reco::Vertex muonVertex = reco::Vertex(transientMuonVertex); + + if (!transientMuonVertex.isValid()) + continue; + std::tuple distanceTuple = + getDistanceBetweenMuonTracks(muonTrack_i, muonTrack_j, magneticField); + // if dca status is good but dca is more than 15 cm + if (std::get<1>(distanceTuple) && std::get<0>(distanceTuple) > 15) + continue; + + nDSADSAVertices++; + + vertexIsValid["DSADSA"].push_back(transientMuonVertex.isValid()); + + std::pair vxy_ = getVxy(muonVertex); + vxy["DSADSA"].push_back(vxy_.first); + vxySigma["DSADSA"].push_back(vxy_.second); + std::pair vxyz_ = getVxyz(muonVertex); + vxyz["DSADSA"].push_back(vxyz_.first); + vxyzSigma["DSADSA"].push_back(vxyz_.second); + chi2["DSADSA"].push_back(muonVertex.chi2()); + ndof["DSADSA"].push_back(muonVertex.ndof()); + normChi2["DSADSA"].push_back(muonVertex.normalizedChi2()); + vx["DSADSA"].push_back(muonVertex.x()); + vy["DSADSA"].push_back(muonVertex.y()); + vz["DSADSA"].push_back(muonVertex.z()); + t["DSADSA"].push_back(muonVertex.t()); + vxErr["DSADSA"].push_back(muonVertex.xError()); + vyErr["DSADSA"].push_back(muonVertex.yError()); + vzErr["DSADSA"].push_back(muonVertex.zError()); + tErr["DSADSA"].push_back(muonVertex.tError()); + dR["DSADSA"].push_back(reco::deltaR(muonTrack_i.eta(), muonTrack_j.eta(), muonTrack_i.phi(), muonTrack_j.phi())); + originalMuonIdx1["DSADSA"].push_back(i); + originalMuonIdx2["DSADSA"].push_back(j); + isDSAMuon1["DSADSA"].push_back(1); + isDSAMuon2["DSADSA"].push_back(1); + + displacedTrackIso03Dimuon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_i, muonVertex, beamSpotInput, &muonTrack_j, 0.3)); + displacedTrackIso04Dimuon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_i, muonVertex, beamSpotInput, &muonTrack_j, 0.4)); + displacedTrackIso03Dimuon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, &muonTrack_i, 0.3)); + displacedTrackIso04Dimuon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, &muonTrack_i, 0.4)); + displacedTrackIso03Muon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_i, muonVertex, beamSpotInput, nullptr, 0.3)); + displacedTrackIso04Muon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_i, muonVertex, beamSpotInput, nullptr, 0.4)); + displacedTrackIso03Muon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, nullptr, 0.3)); + displacedTrackIso04Muon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, muonTrack_j, muonVertex, beamSpotInput, nullptr, 0.4)); + + proxDeltaR["DSADSA"].push_back(getProximityDeltaR(muonTrack_i, muonTrack_j, magneticField, propagator)); + + DCA["DSADSA"].push_back(std::get<0>(distanceTuple)); + DCAstatus["DSADSA"].push_back(std::get<1>(distanceTuple)); + DCAx["DSADSA"].push_back(std::get<2>(distanceTuple).x()); + DCAy["DSADSA"].push_back(std::get<2>(distanceTuple).y()); + DCAz["DSADSA"].push_back(std::get<2>(distanceTuple).z()); + + hitsInFrontOfVert1["DSADSA"].push_back(-1); + missHitsAfterVert1["DSADSA"].push_back(-1); + hitsInFrontOfVert2["DSADSA"].push_back(-1); + missHitsAfterVert2["DSADSA"].push_back(-1); + + reco::TransientTrack refittedTrack_i = transientMuonVertex.refittedTrack(muonTransientTrack_i); + reco::TransientTrack refittedTrack_j = transientMuonVertex.refittedTrack(muonTransientTrack_j); + std::vector refittedTracks = {refittedTrack_i, refittedTrack_j}; + for (size_t k = 0; k < refittedTracks.size(); ++k) { + const auto& refittedTrack = refittedTracks[k]; + refittedTrackIsDSAMuon["DSADSA"].push_back(0); + refittedTrackOriginalIdx["DSADSA"].push_back(i); + refittedTrackPt["DSADSA"].push_back(refittedTrack.track().pt()); + refittedTrackPtErr["DSADSA"].push_back(refittedTrack.track().ptError()); + refittedTrackPx["DSADSA"].push_back(refittedTrack.track().px()); + refittedTrackPy["DSADSA"].push_back(refittedTrack.track().py()); + refittedTrackPz["DSADSA"].push_back(refittedTrack.track().pz()); + refittedTrackEta["DSADSA"].push_back(refittedTrack.track().eta()); + refittedTrackEtaErr["DSADSA"].push_back(refittedTrack.track().etaError()); + refittedTrackPhi["DSADSA"].push_back(refittedTrack.track().phi()); + refittedTrackPhiErr["DSADSA"].push_back(refittedTrack.track().phiError()); + refittedTrackCharge["DSADSA"].push_back(refittedTrack.track().charge()); + refittedTrackNormChi2["DSADSA"].push_back(refittedTrack.normalizedChi2()); + refittedTrackNdof["DSADSA"].push_back(refittedTrack.ndof()); + refittedTrackChi2["DSADSA"].push_back(refittedTrack.chi2()); + + refittedTrackDzPV["DSADSA"].push_back(refittedTrack.track().dz(pv.position())); + refittedTrackDzPVErr["DSADSA"].push_back(std::hypot(refittedTrack.track().dzError(), pv.zError())); + TrajectoryStateClosestToPoint trajectoryPV_i = refittedTrack.trajectoryStateClosestToPoint(primaryVertex); + refittedTrackDxyPVTraj["DSADSA"].push_back(trajectoryPV_i.perigeeParameters().transverseImpactParameter()); + refittedTrackDxyPVTrajErr["DSADSA"].push_back(trajectoryPV_i.perigeeError().transverseImpactParameterError()); + GlobalVector muonRefTrackDir_i( + refittedTrack.track().px(), refittedTrack.track().py(), refittedTrack.track().pz()); + refittedTrackDxyPVSigned["DSADSA"].push_back( + IPTools::signedTransverseImpactParameter(refittedTrack, muonRefTrackDir_i, pv).second.value()); + refittedTrackDxyPVSignedErr["DSADSA"].push_back( + IPTools::signedTransverseImpactParameter(refittedTrack, muonRefTrackDir_i, pv).second.error()); + refittedTrackIp3DPVSigned["DSADSA"].push_back( + IPTools::signedImpactParameter3D(refittedTrack, muonRefTrackDir_i, pv).second.value()); + refittedTrackIp3DPVSignedErr["DSADSA"].push_back( + IPTools::signedImpactParameter3D(refittedTrack, muonRefTrackDir_i, pv).second.error()); + + refittedTrackIdx["DSADSA"].push_back(ddRefittedTrackIdx_counter); + if (k == 0) + refittedTrackIdx1["DSADSA"].push_back(ddRefittedTrackIdx_counter); + if (k == 1) + refittedTrackIdx2["DSADSA"].push_back(ddRefittedTrackIdx_counter); + ddRefittedTrackIdx_counter++; + } + + refittedTrackIso03Dimuon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, &refittedTrack_j.track(), 0.3)); + refittedTrackIso04Dimuon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, &refittedTrack_j.track(), 0.4)); + refittedTrackIso03Dimuon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, &refittedTrack_i.track(), 0.3)); + refittedTrackIso04Dimuon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, &refittedTrack_i.track(), 0.4)); + refittedTrackIso03Muon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, nullptr, 0.3)); + refittedTrackIso04Muon1["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_i.track(), muonVertex, beamSpotInput, nullptr, 0.4)); + refittedTrackIso03Muon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, nullptr, 0.3)); + refittedTrackIso04Muon2["DSADSA"].push_back(getDisplacedTrackerIsolation( + generalTracks, refittedTrack_j.track(), muonVertex, beamSpotInput, nullptr, 0.4)); + } + } + + auto patVertexTab = std::make_unique(nPatPatVertices, "PatMuonVertex", false, false); + auto patdsaVertexTab = std::make_unique(nPatDSAVertices, "PatDSAMuonVertex", false, false); + auto dsaVertexTab = std::make_unique(nDSADSAVertices, "DSAMuonVertex", false, false); + + std::map> vertexTables; + vertexTables["PATPAT"] = std::move(patVertexTab); + vertexTables["PATDSA"] = std::move(patdsaVertexTab); + vertexTables["DSADSA"] = std::move(dsaVertexTab); + + for (const auto& [key, table] : vertexTables) { + table->addColumn("isValid", vertexIsValid[key], ""); + table->addColumn("vxy", vxy[key], ""); + table->addColumn("vxySigma", vxySigma[key], ""); + table->addColumn("vxyz", vxyz[key], ""); + table->addColumn("vxyzSigma", vxyzSigma[key], ""); + table->addColumn("chi2", chi2[key], ""); + table->addColumn("ndof", ndof[key], ""); + table->addColumn("normChi2", normChi2[key], ""); + table->addColumn("vx", vx[key], ""); + table->addColumn("vy", vy[key], ""); + table->addColumn("vz", vz[key], ""); + table->addColumn("t", t[key], ""); + table->addColumn("vxErr", vxErr[key], ""); + table->addColumn("vyErr", vyErr[key], ""); + table->addColumn("vzErr", vzErr[key], ""); + table->addColumn("tErr", tErr[key], ""); + table->addColumn("dR", dR[key], ""); + table->addColumn("originalMuonIdx1", originalMuonIdx1[key], ""); + table->addColumn("originalMuonIdx2", originalMuonIdx2[key], ""); + table->addColumn("isDSAMuon1", isDSAMuon1[key], ""); + table->addColumn("isDSAMuon2", isDSAMuon2[key], ""); + table->addColumn("displacedTrackIso03Dimuon1", displacedTrackIso03Dimuon1[key], ""); + table->addColumn("displacedTrackIso04Dimuon1", displacedTrackIso04Dimuon1[key], ""); + table->addColumn("displacedTrackIso03Dimuon2", displacedTrackIso03Dimuon2[key], ""); + table->addColumn("displacedTrackIso04Dimuon2", displacedTrackIso04Dimuon2[key], ""); + table->addColumn("displacedTrackIso03Muon1", displacedTrackIso03Muon1[key], ""); + table->addColumn("displacedTrackIso04Muon1", displacedTrackIso04Muon1[key], ""); + table->addColumn("displacedTrackIso03Muon2", displacedTrackIso03Muon2[key], ""); + table->addColumn("displacedTrackIso04Muon2", displacedTrackIso04Muon2[key], ""); + table->addColumn("dRprox", proxDeltaR[key], ""); + table->addColumn("dca", DCA[key], ""); + table->addColumn("dcaStatus", DCAstatus[key], ""); + table->addColumn("dcax", DCAx[key], ""); + table->addColumn("dcay", DCAy[key], ""); + table->addColumn("dcaz", DCAz[key], ""); + table->addColumn("hitsInFrontOfVert1", hitsInFrontOfVert1[key], ""); + table->addColumn("hitsInFrontOfVert2", hitsInFrontOfVert2[key], ""); + table->addColumn("missHitsAfterVert1", missHitsAfterVert1[key], ""); + table->addColumn("missHitsAfterVert2", missHitsAfterVert2[key], ""); + table->addColumn("refittedTrackIdx1", refittedTrackIdx1[key], ""); + table->addColumn("refittedTrackIdx2", refittedTrackIdx2[key], ""); + table->addColumn("refittedTrackIso03Dimuon1", refittedTrackIso03Dimuon1[key], ""); + table->addColumn("refittedTrackIso04Dimuon1", refittedTrackIso04Dimuon1[key], ""); + table->addColumn("refittedTrackIso03Dimuon2", refittedTrackIso03Dimuon2[key], ""); + table->addColumn("refittedTrackIso04Dimuon2", refittedTrackIso04Dimuon2[key], ""); + table->addColumn("refittedTrackIso03Muon1", refittedTrackIso03Muon1[key], ""); + table->addColumn("refittedTrackIso04Muon1", refittedTrackIso04Muon1[key], ""); + table->addColumn("refittedTrackIso03Muon2", refittedTrackIso03Muon2[key], ""); + table->addColumn("refittedTrackIso04Muon2", refittedTrackIso04Muon2[key], ""); + } + + iEvent.put(std::move(vertexTables["PATPAT"]), "PatMuonVertex"); + iEvent.put(std::move(vertexTables["PATDSA"]), "PatDSAMuonVertex"); + iEvent.put(std::move(vertexTables["DSADSA"]), "DSAMuonVertex"); + + auto patRefittedTracksTab = + std::make_unique(nPatPatVertices * 2, "PatMuonVertexRefittedTracks", false, false); + auto patdsaRefittedTracksTab = + std::make_unique(nPatDSAVertices * 2, "PatDSAMuonVertexRefittedTracks", false, false); + auto dsaRefittedTracksTab = + std::make_unique(nDSADSAVertices * 2, "DSAMuonVertexRefittedTracks", false, false); + + std::map> refittedTracksTables; + refittedTracksTables["PATPAT"] = std::move(patRefittedTracksTab); + refittedTracksTables["PATDSA"] = std::move(patdsaRefittedTracksTab); + refittedTracksTables["DSADSA"] = std::move(dsaRefittedTracksTab); + + for (const auto& [key, table] : refittedTracksTables) { + table->addColumn("idx", refittedTrackIdx[key], ""); + table->addColumn("isDSAMuon", refittedTrackIsDSAMuon[key], ""); + table->addColumn("originalMuonIdx", refittedTrackOriginalIdx[key], ""); + table->addColumn("pt", refittedTrackPt[key], ""); + table->addColumn("ptErr", refittedTrackPtErr[key], ""); + table->addColumn("px", refittedTrackPx[key], ""); + table->addColumn("py", refittedTrackPy[key], ""); + table->addColumn("pz", refittedTrackPz[key], ""); + table->addColumn("eta", refittedTrackEta[key], ""); + table->addColumn("etaErr", refittedTrackEtaErr[key], ""); + table->addColumn("phi", refittedTrackPhi[key], ""); + table->addColumn("phiErr", refittedTrackPhiErr[key], ""); + table->addColumn("charge", refittedTrackCharge[key], ""); + table->addColumn("normChi2", refittedTrackNormChi2[key], ""); + table->addColumn("ndof", refittedTrackNdof[key], ""); + table->addColumn("chi2", refittedTrackChi2[key], ""); + table->addColumn("dzPV", refittedTrackDzPV[key], ""); + table->addColumn("dzPVErr", refittedTrackDzPVErr[key], ""); + table->addColumn("dxyPVTraj", refittedTrackDxyPVTraj[key], ""); + table->addColumn("dxyPVTrajErr", refittedTrackDxyPVTrajErr[key], ""); + table->addColumn("dxyPVSigned", refittedTrackDxyPVSigned[key], ""); + table->addColumn("dxyPVSignedErr", refittedTrackDxyPVSignedErr[key], ""); + table->addColumn("ip3DPVSigned", refittedTrackIp3DPVSigned[key], ""); + table->addColumn("ip3DPVSignedErr", refittedTrackIp3DPVSignedErr[key], ""); + } + + iEvent.put(std::move(refittedTracksTables["PATPAT"]), "PatMuonVertexRefittedTracks"); + iEvent.put(std::move(refittedTracksTables["PATDSA"]), "PatDSAMuonVertexRefittedTracks"); + iEvent.put(std::move(refittedTracksTables["DSADSA"]), "DSAMuonVertexRefittedTracks"); +} + +std::pair MuonVertexTableProducer::getVxy(const reco::Vertex muonVertex) const { + float vxy = sqrt(muonVertex.x() * muonVertex.x() + muonVertex.y() * muonVertex.y()); + float vxySigma = (1 / vxy) * sqrt(muonVertex.x() * muonVertex.x() * muonVertex.xError() * muonVertex.xError() + + muonVertex.y() * muonVertex.y() * muonVertex.yError() * muonVertex.yError()); + return std::make_pair(vxy, vxySigma); +} + +std::pair MuonVertexTableProducer::getVxyz(const reco::Vertex muonVertex) const { + float vxyz = + sqrt(muonVertex.x() * muonVertex.x() + muonVertex.y() * muonVertex.y() + muonVertex.z() * muonVertex.z()); + float vxyzSigma = (1 / vxyz) * sqrt(muonVertex.x() * muonVertex.x() * muonVertex.xError() * muonVertex.xError() + + muonVertex.y() * muonVertex.y() * muonVertex.yError() * muonVertex.yError() + + muonVertex.z() * muonVertex.z() * muonVertex.zError() * muonVertex.zError()); + return std::make_pair(vxyz, vxyzSigma); +} + +template +float MuonVertexTableProducer::getDisplacedTrackerIsolation(const std::vector& generalTracks, + const MuonType1& muon_1, + const reco::Vertex muonVertex, + const reco::BeamSpot& beamspot, + const MuonType2* muon_2, + float maxDR, + float minDR, + float maxDz, + float maxDxy) const { + float trackPtSum = 0; + + int nGeneralTracks = generalTracks.size(); + float muonTrack2_pt = 0; + float muonTrack2_minDR = 9999; + + for (int i = 0; i < nGeneralTracks; i++) { + const reco::Track& generalTrack = (generalTracks)[i]; + + // Muon POG Tracker Isolation recommendation + float dR = deltaR(muon_1.eta(), muon_1.phi(), generalTrack.eta(), generalTrack.phi()); + if (dR > maxDR) + continue; + if (abs(generalTrack.vz() - muonVertex.z()) > maxDz) + continue; + if (generalTrack.dxy(beamspot) > maxDxy) + continue; + if (dR < minDR) + continue; + + // Determine if track belongs to other muon and get pt of the track + // Only if muon is given as input + if (muon_2 != nullptr) { + float dR_2 = deltaR(muon_2->eta(), muon_2->phi(), generalTrack.eta(), generalTrack.phi()); + if (dR_2 < minDR && dR_2 < muonTrack2_minDR) { + muonTrack2_pt = generalTrack.pt(); + muonTrack2_minDR = dR_2; + } + } + + trackPtSum += generalTrack.pt(); + } + + // Remove pt of track that belongs to other muon + trackPtSum -= muonTrack2_pt; + + float ptRatio = trackPtSum / muon_1.pt(); + return ptRatio; +} + +/** +* Proximity match based on EXO-23-010 +* Calculating deltaR between the innermost hit of the DSAMuon trackRef +* and the extracted closest position of PATMuon track +**/ +template +float MuonVertexTableProducer::getProximityDeltaR(const MuonType& track, + const MuonType& trackRef, + const edm::ESHandle& magneticField, + const edm::ESHandle& propagator) const { + FreeTrajectoryState trajState(GlobalPoint(track.vx(), track.vy(), track.vz()), + GlobalVector(track.px(), track.py(), track.pz()), + track.charge(), + magneticField.product()); + + GlobalPoint refPos(trackRef.innerPosition().x(), trackRef.innerPosition().y(), trackRef.innerPosition().z()); + FreeTrajectoryState trajStatePCA(propagator->propagate(trajState, refPos)); + + float dR = deltaR(trajStatePCA.position().eta(), + trajStatePCA.position().phi(), + trackRef.innerPosition().eta(), + trackRef.innerPosition().phi()); + return dR; +} + +/** +* Proximity between the muons based on EXO-23-010 +* Getting Distance of Closest Approach between muon tracks using TwoTrackMinimumDistance +* Returns tuple of distance (float), error of distance (float) and crossing point (GlobalPoint) +**/ +std::tuple MuonVertexTableProducer::getDistanceBetweenMuonTracks( + const reco::Track& track1, const reco::Track& track2, const edm::ESHandle& magneticField) const { + TwoTrackMinimumDistance ttmd; + FreeTrajectoryState fts1(GlobalPoint(track1.vx(), track1.vy(), track1.vz()), + GlobalVector(track1.px(), track1.py(), track1.pz()), + track1.charge(), + magneticField.product()); + FreeTrajectoryState fts2(GlobalPoint(track2.vx(), track2.vy(), track2.vz()), + GlobalVector(track2.px(), track2.py(), track2.pz()), + track2.charge(), + magneticField.product()); + bool status = ttmd.calculate(fts1, fts2); + if (!status) + return std::tuple(-999.f, status, GlobalPoint(-999.f, -999.f, -999.f)); + return std::make_tuple(ttmd.distance(), status, ttmd.crossingPoint()); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(MuonVertexTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc b/PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc new file mode 100644 index 0000000000000..014b5fb433d0f --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc @@ -0,0 +1,390 @@ +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/NanoAOD/interface/FlatTable.h" +#include "DataFormats/MuonDetId/interface/CSCTriggerNumbering.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/Math/interface/PtEtaPhiMass.h" +#include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h" +#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h" +#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h" +#include "DataFormats/Math/interface/deltaR.h" + +#include "Geometry/CSCGeometry/interface/CSCGeometry.h" +#include "Geometry/DTGeometry/interface/DTGeometry.h" +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" + +#include "fastjet/ClusterSequence.hh" +#include "fastjet/JetDefinition.hh" +#include "fastjet/PseudoJet.hh" + +using RecHitRef = edm::Ref; +using RecHitRefVector = edm::RefVector; + +class cscMDSshowerTableProducer : public edm::global::EDProducer<> { +public: + cscMDSshowerTableProducer(const edm::ParameterSet& iConfig) + : geometryToken_(esConsumes()), + dtgeometryToken_(esConsumes()), + rpcgeometryToken_(esConsumes()), + inputToken_(consumes(iConfig.getParameter("recHitLabel"))), + dtSegmentToken_(consumes(iConfig.getParameter("segmentLabel"))), + rpchitToken_(consumes(iConfig.getParameter("rpcLabel"))), + rParam_(iConfig.getParameter("rParam")), + nRechitMin_(iConfig.getParameter("nRechitMin")), + nStationThres_(iConfig.getParameter("nStationThres")), + stripErr_(iConfig.getParameter("stripErr")), + wireError_(iConfig.getParameter("wireError")), + pruneCut_(iConfig.getParameter("pruneCut")), + name_(iConfig.getParameter("name")) { + produces(name_ + "Rechits"); + produces(name_); + } + + ~cscMDSshowerTableProducer() override {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("recHitLabel")->setComment("input cscRechit collection"); + desc.add("segmentLabel")->setComment("input dt segment collection for veto"); + desc.add("rpcLabel")->setComment("input rpcRechit collection for veto"); + desc.add("rParam", 0.4); + desc.add("nRechitMin", 50); + desc.add("nStationThres", 10); + desc.add("stripErr", 7.0); + desc.add("wireError", 8.6); + desc.add("pruneCut", 9.0); + desc.add("name", "cscRechits")->setComment("name of the output collection"); + descriptions.add("cscMDSshowerTable", desc); + } + float getWeightedTime(RecHitRefVector rechits) const; + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + const edm::ESGetToken geometryToken_; + const edm::ESGetToken dtgeometryToken_; + const edm::ESGetToken rpcgeometryToken_; + edm::EDGetTokenT inputToken_; + edm::EDGetTokenT dtSegmentToken_; + edm::EDGetTokenT rpchitToken_; + const double rParam_; + const int nRechitMin_; // min number of rechits + const int nStationThres_; // min number of rechits to count towards nStation + const double stripErr_, wireError_, pruneCut_; //constants for CSC time + const std::string name_; +}; + +//From: https://github.com/cms-sw/cmssw/blob/master/RecoMuon/MuonIdentification/src/CSCTimingExtractor.cc#L165-L234 +float cscMDSshowerTableProducer::getWeightedTime(RecHitRefVector rechits) const { + bool modified = false; + double totalWeightTimeVtx = 0; + float timeVtx = 0; + + do { + modified = false; + totalWeightTimeVtx = 0; + timeVtx = 0; + for (auto const& rechit : rechits) { + timeVtx += rechit->wireTime() * 1. / (wireError_ * wireError_); + timeVtx += rechit->tpeak() * 1. / (stripErr_ * stripErr_); + totalWeightTimeVtx += 1. / (wireError_ * wireError_); + totalWeightTimeVtx += 1. / (stripErr_ * stripErr_); + } + timeVtx /= totalWeightTimeVtx; + + // cut away outliers + double diff_tvtx_strip, diff_tvtx_wire; + double chimax = 0.0; + int tmmax = 0; + for (size_t i = 0; i < rechits.size(); ++i) { + const auto& rechit = rechits[i]; + //diff_tvtx = (cscHits[i].time - timeVtx) * (cscHits[i].time - timeVtx) * cscHits[i].error; + diff_tvtx_strip = (rechit->tpeak() - timeVtx) * (rechit->tpeak() - timeVtx) * 1. / (stripErr_ * stripErr_); + diff_tvtx_wire = (rechit->wireTime() - timeVtx) * (rechit->wireTime() - timeVtx) * 1. / (wireError_ * wireError_); + + if ((diff_tvtx_strip > chimax) || (diff_tvtx_wire > chimax)) { + tmmax = i; + chimax = std::max(diff_tvtx_strip, diff_tvtx_wire); + } + } + // cut away the outliers and repeat + if (chimax > pruneCut_) { + rechits.erase(rechits.begin() + tmmax); + modified = true; + } + } while (modified); + + return timeVtx; +} + +void cscMDSshowerTableProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + auto const& geo = iSetup.getData(geometryToken_); + auto const& dt_geo = iSetup.getData(dtgeometryToken_); + auto const& rpc_geo = iSetup.getData(rpcgeometryToken_); + auto const& rechits = iEvent.get(inputToken_); + auto const& segments = iEvent.get(dtSegmentToken_); + auto const& rpchits = iEvent.get(rpchitToken_); + + std::set unique_ids; + std::vector fjInput; + + edm::RefVector inputs; + + fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, rParam_); + + int recIt = 0; + for (auto const& rechit : rechits) { + LocalPoint recHitLocalPosition = rechit.localPosition(); + auto detid = rechit.cscDetId(); + auto thischamber = geo.chamber(detid); + if (thischamber) { + GlobalPoint globalPosition = thischamber->toGlobal(recHitLocalPosition); + float x = globalPosition.x(); + float y = globalPosition.y(); + float z = globalPosition.z(); + RecHitRef ref = RecHitRef(&rechits, recIt); + inputs.push_back(ref); + fjInput.push_back(fastjet::PseudoJet(x, y, z, globalPosition.mag())); + fjInput.back().set_user_index(recIt); + } + recIt++; + } + fastjet::ClusterSequence clus_seq(fjInput, jet_def); + + //keep all the clusters + double ptmin = 0.0; + std::vector fjJets = clus_seq.inclusive_jets(ptmin); + + // Constituent rechit fields + std::vector cscRechitsX, cscRechitsY, cscRechitsZ, cscRechitsPhi, cscRechitsEta, cscRechitsE, cscRechitsTpeak, + cscRechitsTwire; + std::vector cscRechitsNStrips, cscRechitsHitWire, cscRechitsWGroupsBX, cscRechitsNWireGroups, cscRechitsQuality, + cscRechitsChamber, cscRechitsIChamber, cscRechitsStation; + ; + + // MDS fields + std::vector clsX, clsY, clsZ, clsPhi, clsEta, clsTime, clsTimeSpread, clsTimeWeighted, clsTimeSpreadWeighted, + clsAvgStation; + std::vector clsSize, clsNstation, cls_nME11, cls_nME12, clsUniqueChamber, cls_nMB1dtSeg, cls_nRE12hit, + cls_nRB1hit; + + for (auto const& fjJet : fjJets) { + // skip if the cluster has too few rechits + if (int(fjJet.constituents().size()) < nRechitMin_) + continue; + // get the constituents from fastjet + RecHitRefVector rechits; + for (auto const& constituent : fjJet.constituents()) { + auto index = constituent.user_index(); + if (index >= 0 && static_cast(index) < inputs.size()) { + rechits.push_back(inputs[index]); + } + } + + //Derive cluster properties + int nME12 = 0, nME11 = 0, nMB1dtSeg = 0, nRE12hit = 0, nRB1hit = 0; + int nStation = 0; + int totStation = 0; + float avgStation = 0.0; + float timeSpread = 0.0; + float time = 0.0; + float time_strip = 0.0; // for timeSpread calculation + double timeWeighted = 0.0; + float timeSpreadWeighted = 0.0; + + std::map station_count_map; + + //fill rechits fields + for (auto const& rechit : rechits) { + LocalPoint recHitLocalPosition = rechit->localPosition(); + auto detid = rechit->cscDetId(); + unique_ids.insert(detid.chamberId()); + auto thischamber = geo.chamber(detid); + int endcap = CSCDetId::endcap(detid) == 1 ? 1 : -1; + if (thischamber) { + GlobalPoint globalPosition = thischamber->toGlobal(recHitLocalPosition); + + cscRechitsX.push_back(globalPosition.x()); + cscRechitsY.push_back(globalPosition.y()); + cscRechitsZ.push_back(globalPosition.z()); + cscRechitsPhi.push_back(globalPosition.phi()); + cscRechitsEta.push_back(globalPosition.eta()); + cscRechitsE.push_back(rechit->energyDepositedInLayer()); //not saved + cscRechitsTpeak.push_back(rechit->tpeak()); + cscRechitsTwire.push_back(rechit->wireTime()); + cscRechitsQuality.push_back(rechit->quality()); + + int stationRing = (CSCDetId::station(detid) * 10 + CSCDetId::ring(detid)); + if (CSCDetId::ring(detid) == 4) + stationRing = (CSCDetId::station(detid) * 10 + 1); // ME1/a has ring==4 + + cscRechitsChamber.push_back(endcap * stationRing); + cscRechitsIChamber.push_back(CSCDetId::chamber(detid)); + cscRechitsStation.push_back(endcap * CSCDetId::station(detid)); + cscRechitsNStrips.push_back(rechit->nStrips()); + cscRechitsHitWire.push_back(rechit->hitWire()); + cscRechitsWGroupsBX.push_back(rechit->wgroupsBX()); + cscRechitsNWireGroups.push_back(rechit->nWireGroups()); + + //compute for cluster fields + station_count_map[CSCDetId::station(detid)]++; + if (stationRing == 11) + nME11++; + if (stationRing == 12) + nME12++; + time += (rechit->tpeak() + rechit->wireTime()); + time_strip += rechit->tpeak(); + } + } + //station statistics + std::map::iterator it; + for (auto const& [station, count] : station_count_map) { + if (count >= nStationThres_) { + nStation++; + avgStation += station * count; + totStation += count; + } + } + if (totStation != 0) { + avgStation = avgStation / totStation; + } + float invN = 1.f / rechits.size(); + time = (time / 2.f) * invN; + time_strip = time_strip * invN; + + // https://github.com/cms-sw/cmssw/blob/master/RecoMuon/MuonIdentification/src/CSCTimingExtractor.cc#L165-L234 + for (auto& rechit : rechits) { + timeSpread += (rechit->tpeak() - time_strip) * (rechit->tpeak() - time_strip); + } + timeSpread = std::sqrt(timeSpread * invN); + + float i_clsEta = etaFromXYZ(fjJet.px() * invN, fjJet.py() * invN, fjJet.pz() * invN); + float i_clsPhi = std::atan2(fjJet.py() * invN, fjJet.px() * invN); + + //MB1 DT seg + for (auto const& segment : segments) { + LocalPoint localPosition = segment.localPosition(); + auto geoid = segment.geographicalId(); + DTChamberId dtdetid = DTChamberId(geoid); + auto thischamber = dt_geo.chamber(dtdetid); + if (thischamber) { + GlobalPoint globalPosition = thischamber->toGlobal(localPosition); + float eta = globalPosition.eta(); + float phi = globalPosition.phi(); + if (dtdetid.station() == 1 && reco::deltaR(eta, phi, i_clsEta, i_clsPhi) < 0.4) + nMB1dtSeg++; + } + } + //RPC hits + for (auto const& rechit : rpchits) { + LocalPoint recHitLocalPosition = rechit.localPosition(); + auto geoid = rechit.geographicalId(); + RPCDetId rpcdetid = RPCDetId(geoid); + auto thischamber = rpc_geo.chamber(rpcdetid); + + if (thischamber) { + GlobalPoint globalPosition = thischamber->toGlobal(recHitLocalPosition); + float eta = globalPosition.eta(); + float phi = globalPosition.phi(); + + //RE12 hits + if (rpcdetid.station() == 1 && rpcdetid.ring() == 2 && abs(rpcdetid.region()) == 1 && + reco::deltaR(eta, phi, i_clsEta, i_clsPhi) < 0.4) + nRE12hit++; + //RB1 hits + if (rpcdetid.station() == 1 && rpcdetid.region() == 0 && reco::deltaR(eta, phi, i_clsEta, i_clsPhi) < 0.4) + nRB1hit++; + } + } + + //fill cluster fields + clsSize.push_back(rechits.size()); + // cluster position is the average position of the constituent rechits + clsX.push_back(fjJet.px() * invN); + clsY.push_back(fjJet.py() * invN); + clsZ.push_back(fjJet.pz() * invN); + clsEta.push_back(i_clsEta); + clsPhi.push_back(i_clsPhi); + + clsTime.push_back(time); + clsTimeSpread.push_back(timeSpread); + cls_nME11.push_back(nME11); + cls_nME12.push_back(nME12); + clsNstation.push_back(nStation); + clsAvgStation.push_back(avgStation); + clsUniqueChamber.push_back(unique_ids.size()); + cls_nMB1dtSeg.push_back(nMB1dtSeg); + cls_nRE12hit.push_back(nRE12hit); + cls_nRB1hit.push_back(nRB1hit); + + //cluster time weighted with unc. and pruned outliers + timeWeighted = cscMDSshowerTableProducer::getWeightedTime(rechits); + + for (auto& rechit : rechits) { + timeSpreadWeighted += (timeWeighted - rechit->tpeak()) * (timeWeighted - rechit->tpeak()); + } + timeSpreadWeighted = std::sqrt(timeSpreadWeighted * invN); + + clsTimeWeighted.push_back(timeWeighted); + clsTimeSpreadWeighted.push_back(timeSpreadWeighted); + } + auto cscRechitTab = std::make_unique(cscRechitsX.size(), name_ + "Rechits", false, false); + + cscRechitTab->addColumn("X", cscRechitsX, "csc rechit X"); + cscRechitTab->addColumn("Y", cscRechitsY, "csc rechit Y"); + cscRechitTab->addColumn("Z", cscRechitsZ, "csc rechit Z"); + cscRechitTab->addColumn("Phi", cscRechitsPhi, "csc rechit Phi"); + cscRechitTab->addColumn("Eta", cscRechitsEta, "csc rechit Eta"); + cscRechitTab->addColumn("E", cscRechitsE, "csc rechit Energy deposited in layer"); + cscRechitTab->addColumn("Tpeak", cscRechitsTpeak, "csc rechit time from cathode"); + cscRechitTab->addColumn("Twire", cscRechitsTwire, "csc rechit time from anode"); + cscRechitTab->addColumn("Quality", cscRechitsQuality, "csc rechit quality"); + cscRechitTab->addColumn("Chamber", cscRechitsChamber, "csc rechit station-Ring"); + cscRechitTab->addColumn("IChamber", cscRechitsIChamber, "csc rechit chamber in ring"); + cscRechitTab->addColumn("Station", cscRechitsStation, "csc rechit station"); + cscRechitTab->addColumn("NStrips", cscRechitsNStrips, "csc rechit nstrips"); + cscRechitTab->addColumn("WGroupsBX", cscRechitsWGroupsBX, "csc rechit wire group BX"); + cscRechitTab->addColumn("HitWire", cscRechitsHitWire, "csc rechit hit wire"); + cscRechitTab->addColumn("NWireGroups", cscRechitsNWireGroups, "csc rechit n wire groups"); + + iEvent.put(std::move(cscRechitTab), name_ + "Rechits"); + + auto clsTab = std::make_unique(clsSize.size(), name_, false, false); + + clsTab->addColumn("size", clsSize, "cluster Size"); + clsTab->addColumn("x", clsX, "cluster X"); + clsTab->addColumn("y", clsY, "cluster Y"); + clsTab->addColumn("z", clsZ, "cluster Z"); + clsTab->addColumn("phi", clsPhi, "cluster Phi"); + clsTab->addColumn("eta", clsEta, "cluster Eta"); + clsTab->addColumn("time", clsTime, "cluster Time"); + clsTab->addColumn("timeSpread", clsTimeSpread, "cluster TimeSpread"); + clsTab->addColumn("timeWeighted", clsTimeWeighted, "cluster TimeWeighted"); + clsTab->addColumn("timeSpreadWeighted", clsTimeSpreadWeighted, "cluster TimeSpreadWeighted"); + clsTab->addColumn("nStation", clsNstation, "cluster nStation"); + clsTab->addColumn("uniqueChamber", clsUniqueChamber, "cluster unique chambers"); + clsTab->addColumn("avgStation", clsAvgStation, "cluster AvgStation"); + clsTab->addColumn("nME11", cls_nME11, "cluster nME11"); + clsTab->addColumn("nME12", cls_nME12, "cluster nME12"); + clsTab->addColumn("nMB1dtSeg", cls_nMB1dtSeg, "cluster nMB1dtSeg"); + clsTab->addColumn("nRE12hit", cls_nRE12hit, "cluster nRE12hit"); + clsTab->addColumn("nRB1hit", cls_nRB1hit, "cluster nRB1hit"); + + iEvent.put(std::move(clsTab), name_); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(cscMDSshowerTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/dtMDSshowerTableProducer.cc b/PhysicsTools/NanoAOD/plugins/dtMDSshowerTableProducer.cc new file mode 100644 index 0000000000000..28a3a0a1d6eeb --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/dtMDSshowerTableProducer.cc @@ -0,0 +1,339 @@ +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/global/EDProducer.h" + +#include "FWCore/Framework/interface/Event.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/NanoAOD/interface/FlatTable.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/Math/interface/PtEtaPhiMass.h" +#include "DataFormats/DTRecHit/interface/DTRecHitCollection.h" +#include "DataFormats/RPCRecHit/interface/RPCRecHitCollection.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +#include "Geometry/DTGeometry/interface/DTGeometry.h" +#include "Geometry/RPCGeometry/interface/RPCGeometry.h" +#include "Geometry/Records/interface/MuonGeometryRecord.h" + +#include "fastjet/ClusterSequence.hh" +#include "fastjet/JetDefinition.hh" +#include "fastjet/PseudoJet.hh" + +using RecHitRef = edm::Ref; +using RecHitRefVector = edm::RefVector; + +class dtMDSshowerTableProducer : public edm::global::EDProducer<> { +public: + dtMDSshowerTableProducer(const edm::ParameterSet& iConfig) + : dtgeometryToken_(esConsumes()), + rpcgeometryToken_(esConsumes()), + inputToken_(consumes(iConfig.getParameter("recHitLabel"))), + rpchitToken_(consumes(iConfig.getParameter("rpcLabel"))), + rParam_(iConfig.getParameter("rParam")), + nRechitMin_(iConfig.getParameter("nRechitMin")), + nStationThres_(iConfig.getParameter("nStationThres")), + name_(iConfig.getParameter("name")) { + produces(name_ + "Rechits"); + produces(name_); + } + + ~dtMDSshowerTableProducer() override {} + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("recHitLabel")->setComment("input dtRechit collection"); + desc.add("rpcLabel")->setComment("input rpcRechit collection for veto"); + desc.add("rParam", 0.4); + desc.add("nRechitMin", 50); + desc.add("nStationThres", 10); + desc.add("name", "dt1DRecHits")->setComment("name of the output collection"); + descriptions.add("dtMDSshowerTable", desc); + } + float getWeightedTime(RecHitRefVector rechits) const; + +private: + void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override; + + const edm::ESGetToken dtgeometryToken_; + const edm::ESGetToken rpcgeometryToken_; + edm::EDGetTokenT inputToken_; + edm::EDGetTokenT rpchitToken_; + const double rParam_; + const int nRechitMin_; // min number of rechits + const int nStationThres_; // min number of rechits to count towards nStation + const std::string name_; +}; + +void dtMDSshowerTableProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const { + auto const& dt_geo = iSetup.getData(dtgeometryToken_); + auto const& rpc_geo = iSetup.getData(rpcgeometryToken_); + auto const& rechits = iEvent.get(inputToken_); + auto const& rpchits = iEvent.get(rpchitToken_); + + std::set unique_ids; + std::vector fjInput; + + RecHitRefVector inputs; + + fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, rParam_); + + int recIt = 0; + for (auto const& rechit : rechits) { + LocalPoint recHitLocalPosition = rechit.localPosition(); + auto geoid = rechit.geographicalId(); + DTChamberId dtdetid = DTChamberId(geoid); + DTLayerId dtlayerid = DTLayerId(geoid); + + auto thischamber = dt_geo.chamber(dtdetid); + auto thislayer = dt_geo.layer(dtlayerid); + + if (thischamber) { + GlobalPoint globalPosition; + if (thislayer) { + globalPosition = thislayer->toGlobal(recHitLocalPosition); + } else { + globalPosition = thischamber->toGlobal(recHitLocalPosition); + } + float x = globalPosition.x(); + float y = globalPosition.y(); + float z = globalPosition.z(); + RecHitRef ref = RecHitRef(&rechits, recIt); + inputs.push_back(ref); + fjInput.push_back(fastjet::PseudoJet(x, y, z, globalPosition.mag())); + fjInput.back().set_user_index(recIt); + } + recIt++; + } + fastjet::ClusterSequence clus_seq(fjInput, jet_def); + + //keep all the clusters + double ptmin = 0.0; + std::vector fjJets = clus_seq.inclusive_jets(ptmin); + + // Constituent rechit fields + std::vector dtRechitsX, dtRechitsY, dtRechitsZ, dtRechitsPhi, dtRechitsEta; + std::vector dtRechitsLayer, dtRechitsSuperLayer, dtRechitsSector, dtRechitsStation, dtRechitsWheel; + + // MDS fields + std::vector clsX, clsY, clsZ, clsPhi, clsEta, clsAvgStation; + std::vector clsSize, clsNstation, clsWheel, clsUniqueChamber, cls_nRB1hit, cls_nRPC, clsBX; + + for (auto const& fjJet : fjJets) { + // skip if the cluster has too few rechits + if (int(fjJet.constituents().size()) < nRechitMin_) + continue; + // get the constituents from fastjet + RecHitRefVector rechits; + for (auto const& constituent : fjJet.constituents()) { + auto index = constituent.user_index(); + if (index >= 0 && static_cast(index) < inputs.size()) { + rechits.push_back(inputs[index]); + } + } + + //Derive cluster properties + int size_z = 0, size_xy = 0, size = 0; + float avg_x_sl2(0.0), avg_y_sl2(0.0), avg_z_sl2(0.0); + float avg_x(0.0), avg_y(0.0), avg_z(0.0); + int nStation = 0; + int totStation = 0; + float avgStation = 0.0; + + std::map station_count_map; + + //fill rechits fields + for (auto const& rechit : rechits) { + LocalPoint recHitLocalPosition = rechit->localPosition(); + auto geoid = rechit->geographicalId(); + + DTChamberId dtdetid = DTChamberId(geoid); + DTLayerId dtlayerid = DTLayerId(geoid); + + unique_ids.insert(dtdetid); + + auto thischamber = dt_geo.chamber(dtdetid); + auto thislayer = dt_geo.layer(dtlayerid); + + if (thischamber) { + if (thislayer) { + GlobalPoint globalPosition = thislayer->toGlobal(recHitLocalPosition); + dtRechitsX.push_back(globalPosition.x()); + dtRechitsY.push_back(globalPosition.y()); + dtRechitsZ.push_back(globalPosition.z()); + dtRechitsPhi.push_back(globalPosition.phi()); + dtRechitsEta.push_back(globalPosition.eta()); + dtRechitsLayer.push_back(dtlayerid.layer()); + dtRechitsSuperLayer.push_back(dtlayerid.superlayer()); + // use xy-coordinates from SL1/SL3 when available + if (dtlayerid.superlayer() == 1 || dtlayerid.superlayer() == 3) { + avg_x += globalPosition.x(); + avg_y += globalPosition.y(); + avg_z += globalPosition.z(); + size_xy++; + } + // use z-coordinates from SL2 when available + else if (dtlayerid.superlayer() == 2) { + avg_x_sl2 += globalPosition.x(); + avg_y_sl2 += globalPosition.y(); + avg_z_sl2 += globalPosition.z(); + size_z++; + } + } else { + GlobalPoint globalPosition = thischamber->toGlobal(recHitLocalPosition); + dtRechitsX.push_back(globalPosition.x()); + dtRechitsY.push_back(globalPosition.y()); + dtRechitsZ.push_back(globalPosition.z()); + dtRechitsPhi.push_back(globalPosition.phi()); + dtRechitsEta.push_back(globalPosition.eta()); + dtRechitsLayer.push_back(0); //default value; + dtRechitsSuperLayer.push_back(0); //default value + avg_x += globalPosition.x(); + avg_y += globalPosition.y(); + avg_z += globalPosition.z(); + } + dtRechitsSector.push_back(dtdetid.sector()); + dtRechitsStation.push_back(dtdetid.station()); + dtRechitsWheel.push_back(dtdetid.wheel()); + size++; + + //compute for cluster fields + station_count_map[dtdetid.station()]++; + } + } + //station statistics + std::map::iterator it; + for (auto const& [station, count] : station_count_map) { + if (count >= nStationThres_) { + nStation++; + avgStation += station * count; + totStation += count; + } + } + if (totStation != 0) { + avgStation = avgStation / totStation; + } + + //for DT correct position, calculate average Z using sl2 and average XY using sl1/3 + if (size_xy > 0 && size_z > 0) { // both SL1/SL3 and SL2 rechits + avg_x = avg_x / size_xy; + avg_y = avg_y / size_xy; + avg_z = avg_z_sl2 / size_z; + } else if (size_xy == 0 && size_z > 0) // only SL2 rechits + { + avg_x = avg_x_sl2 / size_z; + avg_y = avg_y_sl2 / size_z; + avg_z = avg_z_sl2 / size_z; + } else if (size_xy > 0 && size_z == 0) // no SL2 rechits + { + avg_x = avg_x / size_xy; + avg_y = avg_y / size_xy; + avg_z = avg_z / size_xy; + } else // no SL information + { + avg_x = avg_x / size; + avg_y = avg_y / size; + avg_z = avg_z / size; + } + + float i_clsEta = etaFromXYZ(avg_x, avg_y, avg_z); + float i_clsPhi = std::atan2(avg_y, avg_x); + int i_clsWheel = (std::abs(avg_z) < 126.8) ? 0 + : (avg_z > 126.8 && avg_z < 395.4) ? 1 + : (avg_z < -126.8 && avg_z > -395.4) ? -1 + : (avg_z < 0) ? -2 + : 2; + + // RPC hits + std::map bxCounts; // bx of matched RPC hits : counts + int nRB1hit = 0; + for (auto const& rechit : rpchits) { + LocalPoint recHitLocalPosition = rechit.localPosition(); + auto geoid = rechit.geographicalId(); + RPCDetId rpcdetid = RPCDetId(geoid); + auto thischamber = rpc_geo.chamber(rpcdetid); + + // Only matching RB hits + if (rpcdetid.region() != 0) + continue; + if (thischamber) { + GlobalPoint globalPosition = thischamber->toGlobal(recHitLocalPosition); + + //match to RPC hits with dPhi<0.5 and same wheel in DT + if (reco::deltaPhi(globalPosition.phi(), i_clsPhi) < 0.5 && rpcdetid.ring() == i_clsWheel) { + //RB1 hits + if (rpcdetid.station() == 1) + nRB1hit++; + bxCounts[rechit.BunchX()]++; + } + } + } + // find the mode of BX + int modeBX = 0, maxCount = 0, i_cls_nRPC = 0; + for (const auto& [bx, count] : bxCounts) { + i_cls_nRPC += count; + if (count > maxCount) { + modeBX = bx; + maxCount = count; + } + } + + //fill cluster fields + clsSize.push_back(rechits.size()); + // cluster position is the average position of the constituent rechits + clsX.push_back(avg_x); + clsY.push_back(avg_y); + clsZ.push_back(avg_z); + clsEta.push_back(i_clsEta); + clsPhi.push_back(i_clsPhi); + clsBX.push_back(modeBX); + clsNstation.push_back(nStation); + clsAvgStation.push_back(avgStation); + clsWheel.push_back(i_clsWheel); + clsUniqueChamber.push_back(unique_ids.size()); + cls_nRB1hit.push_back(nRB1hit); + cls_nRPC.push_back(i_cls_nRPC); + } + auto dtRechitTab = std::make_unique(dtRechitsX.size(), name_ + "Rechits", false, false); + + dtRechitTab->addColumn("X", dtRechitsX, "dt rechit X"); + dtRechitTab->addColumn("Y", dtRechitsY, "dt rechit Y"); + dtRechitTab->addColumn("Z", dtRechitsZ, "dt rechit Z"); + dtRechitTab->addColumn("Phi", dtRechitsPhi, "dt rechit Phi"); + dtRechitTab->addColumn("Eta", dtRechitsEta, "dt rechit Eta"); + dtRechitTab->addColumn("Layer", dtRechitsLayer, "dt rechit Layer"); + dtRechitTab->addColumn("SuperLayer", dtRechitsSuperLayer, "dt rechit SuperLayer"); + dtRechitTab->addColumn("Sector", dtRechitsSector, "dt rechit sector"); + dtRechitTab->addColumn("Station", dtRechitsStation, "dt rechit station"); + dtRechitTab->addColumn("Wheel", dtRechitsWheel, "dt rechit nstrips"); + + iEvent.put(std::move(dtRechitTab), name_ + "Rechits"); + + auto clsTab = std::make_unique(clsSize.size(), name_, false, false); + + clsTab->addColumn("size", clsSize, "cluster Size"); + clsTab->addColumn("x", clsX, "cluster X"); + clsTab->addColumn("y", clsY, "cluster Y"); + clsTab->addColumn("z", clsZ, "cluster Z"); + clsTab->addColumn("phi", clsPhi, "cluster Phi"); + clsTab->addColumn("eta", clsEta, "cluster Eta"); + clsTab->addColumn("bx", clsBX, "cluster BX"); + clsTab->addColumn("wheel", clsWheel, "cluster wheel"); + clsTab->addColumn("nStation", clsNstation, "cluster nStation"); + clsTab->addColumn("uniqueChamber", clsUniqueChamber, "cluster unique chambers"); + clsTab->addColumn("avgStation", clsAvgStation, "cluster AvgStation"); + clsTab->addColumn("nRPC", cls_nRPC, "cluster nRPC"); + clsTab->addColumn("nRB1hit", cls_nRB1hit, "cluster nRB1hit"); + + iEvent.put(std::move(clsTab), name_); +} + +#include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(dtMDSshowerTableProducer); diff --git a/PhysicsTools/NanoAOD/python/autoNANO.py b/PhysicsTools/NanoAOD/python/autoNANO.py index e39bcb67d91b2..4cb061ee5dff7 100644 --- a/PhysicsTools/NanoAOD/python/autoNANO.py +++ b/PhysicsTools/NanoAOD/python/autoNANO.py @@ -99,6 +99,9 @@ def expandNanoMapping(seqList, mapping, key): # Custom BTV Nano for SF measurements or tagger training 'BTV': {'sequence': '@PHYS', 'customize': '@PHYS+PhysicsTools/NanoAOD/custom_btv_cff.BTVCustomNanoAOD'}, + # Custom EXO Nano for EXO analyses + 'EXO': {'sequence': '@PHYS', + 'customize': '@PHYS+PhysicsTools/NanoAOD/custom_exo_cff.add_exonanoTables'}, # NANOGEN (from LHE/GEN/AOD) 'GEN': {'sequence': 'PhysicsTools/NanoAOD/nanogen_cff.nanogenSequence', 'customize': 'PhysicsTools/NanoAOD/nanogen_cff.customizeNanoGEN'}, diff --git a/PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py b/PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py new file mode 100644 index 0000000000000..7318d933da82a --- /dev/null +++ b/PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py @@ -0,0 +1,157 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * +from PhysicsTools.NanoAOD.genparticles_cff import * +from PhysicsTools.NanoAOD.taus_cff import * +from PhysicsTools.NanoAOD.muons_cff import * +from PhysicsTools.NanoAOD.jetsAK4_CHS_cff import * +from PhysicsTools.NanoAOD.jetsAK4_Puppi_cff import * +import os + +def add_displacedtauCHSTables(process, isMC): + + process.linkedObjectsCHS = cms.EDProducer("PATObjectCrossLinker", + jets=cms.InputTag("finalJets"), + muons=cms.InputTag("finalMuons"), + electrons=cms.InputTag("finalElectrons"), + lowPtElectrons=cms.InputTag("finalLowPtElectrons"), + taus=cms.InputTag("finalTaus"), + boostedTaus=cms.InputTag("finalBoostedTaus"), + photons=cms.InputTag("finalPhotons"), + vertices=cms.InputTag("slimmedSecondaryVertices") + ) + + del process.updatedJetsWithUserData.userFloats.leadTrackPt + del process.updatedJetsWithUserData.userFloats.leptonPtRelv0 + del process.updatedJetsWithUserData.userFloats.leptonPtRelInvv0 + del process.updatedJetsWithUserData.userFloats.leptonDeltaR + del process.updatedJetsWithUserData.userFloats.vtxPt + del process.updatedJetsWithUserData.userFloats.vtxMass + del process.updatedJetsWithUserData.userFloats.vtx3dL + del process.updatedJetsWithUserData.userFloats.vtx3deL + del process.updatedJetsWithUserData.userFloats.ptD + del process.updatedJetsWithUserData.userFloats.qgl + del process.updatedJetsWithUserData.userFloats.puIdNanoDisc + del process.updatedJetsWithUserData.userFloats.muonSubtrRawPt + del process.updatedJetsWithUserData.userFloats.muonSubtrRawEta + del process.updatedJetsWithUserData.userFloats.muonSubtrRawPhi + + del process.updatedJetsWithUserData.userInts.vtxNtrk + del process.updatedJetsWithUserData.userInts.leptonPdgId + del process.updatedJetsWithUserData.userInts.puIdNanoId + + print(process.updatedJetsWithUserData.dumpPython()) + + # + # Customize jetTable + # + process.jetTable.src = cms.InputTag("linkedObjectsCHS","jets") + process.jetTable.name = "JetCHS" # Change collection name from "Jet" ->" JetCHS" + + # + # Remove these tagger branches since for CHS, we just want to store ParticleNet. + # Remove also branches related to object linking. It is only done for AK4 Puppi. + # + for varName in process.jetTable.variables.parameterNames_(): + if "btagDeepFlav" in varName or "btagRobustParT" in varName or "btagUParT" in varName: + delattr(process.jetTable.variables, varName) + if "UParTAK4Reg" in varName: + delattr(process.jetTable.variables, varName) + if "svIdx" in varName or "muonIdx" in varName or "electronIdx" in varName: + delattr(process.jetTable.variables, varName) + if "nSVs" in varName or "nElectrons" in varName or "nMuons" in varName: + delattr(process.jetTable.variables, varName) + + del process.jetTable.variables.muonSubtrFactor + del process.jetTable.variables.muonSubtrDeltaEta + del process.jetTable.variables.muonSubtrDeltaPhi + del process.jetTable.variables.qgl + del process.jetTable.variables.puIdDisc + del process.jetTable.variables.puId + + del process.jetTable.externalVariables.bRegCorr + del process.jetTable.externalVariables.bRegRes + del process.jetTable.externalVariables.cRegCorr + del process.jetTable.externalVariables.cRegRes + + process.jetUserDataTask = cms.Task( + process.jercVars, + ) + process.nanoTableTaskCommon.add(process.jetUserDataTask) + + + + + ## displaced tau part + file = "RecoTauTag/TrainingFiles/DisplacedTauId/particlenet_v1_a27159734e304ea4b7f9e0042baa9e22.pb" + process.options = cms.untracked.PSet( + numberOfThreads = cms.untracked.uint32(4), # Global thread count + numberOfStreams = cms.untracked.uint32(4), # Should match threads + ) + + process.disTauTag = cms.EDProducer( + "DisTauTag", + graphPath = cms.FileInPath(file), + jets = cms.InputTag("linkedObjectsCHS","jets"), + pfCandidates = cms.InputTag('packedPFCandidates'), + save_inputs = cms.bool(False), + batchSize = cms.uint32(8), + #numThreads = cms.untracked.uint32(4) + allowUnscheduled = cms.untracked.bool(True) + ) + + process.jetImpactParameters = cms.EDProducer( + "JetImpactParameters", + jets = cms.InputTag("linkedObjectsCHS","jets"), + pfCandidates = cms.InputTag('packedPFCandidates'), + deltaRMax = cms.double(0.4) + ) + + + + d_disTauTagVars = { + "disTauTag_score0": ExtVar("disTauTag:score0" , float, doc = "Score 0"), + "disTauTag_score1": ExtVar("disTauTag:score1" , float, doc = "Score 1"), + "dxy": ExtVar("jetImpactParameters:jetDxy", float, doc = "leadingPtPFCand_dxy which is within dR=0.4 and charged/hasTrackDetails"), + "dz": ExtVar("jetImpactParameters:jetDz", float, doc = "leadingPtPFCand_dz which is within dR=0.4 and charged/hasTrackDetails"), + "dxyerror": ExtVar("jetImpactParameters:jetDxyError", float, doc = "leadingPtPFCand_dxyerror which is within dR=0.4 and charged/hasTrackDetails"), + "dzerror": ExtVar("jetImpactParameters:jetDzError", float, doc = "leadingPtPFCand_dzerror which is within dR=0.4 and charged/hasTrackDetails"), + "charge": ExtVar("jetImpactParameters:jetCharge", float, doc = "leadingPtPFCand_charge which is within dR=0.4 and charged/hasTrackDetails"), + } + + #if useCHSJets: + process.jetTable.externalVariables = process.jetTable.externalVariables.clone(**d_disTauTagVars) + ## for puppi jets, use this! + #else: + # process.jetPuppiTable.externalVariables = process.jetPuppiTable.externalVariables.clone(**d_disTauTagVars) + + process.jetTask = cms.Task( + process.jetCorrFactorsNano, + process.updatedJets, + process.linkedObjectsCHS, + process.updatedJetsWithUserData, + process.finalJets, + process.disTauTag, + process.jetImpactParameters + + ) + + process.nanoTableTaskCommon.add(process.jetTask) + + process.jetTablesTask = cms.Task( + process.jetTable + ) + process.nanoTableTaskCommon.add(process.jetTablesTask) + + # + # Only for MC + # + if isMC: + process.jetCHSMCTable = process.jetMCTable.clone( + src = process.jetTable.src, + name = process.jetTable.name + ) + process.jetMCTask.add(process.jetCHSMCTable) + + return process + + diff --git a/PhysicsTools/NanoAOD/python/custom_exo_cff.py b/PhysicsTools/NanoAOD/python/custom_exo_cff.py new file mode 100644 index 0000000000000..9cedea34d92a1 --- /dev/null +++ b/PhysicsTools/NanoAOD/python/custom_exo_cff.py @@ -0,0 +1,315 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * + +# displaced tau +from PhysicsTools.NanoAOD.custom_displacedtau_cff import * + +unpackedTracksAndVertices = cms.EDProducer('PATTrackAndVertexUnpacker', + slimmedVertices = cms.InputTag("offlineSlimmedPrimaryVertices"), + slimmedSecondaryVertices = cms.InputTag("slimmedSecondaryVertices"), + additionalTracks= cms.InputTag("lostTracks"), + packedCandidates = cms.InputTag("packedPFCandidates")) + +# ref: https://github.com/cms-sw/cmssw/blob/2bed69b1658e4deeaef914e462741919e9183be3/RecoVertex/AdaptiveVertexFinder/plugins/InclusiveVertexFinder.h#L48 + +displacedInclusiveVertexFinder = cms.EDProducer("InclusiveVertexFinder", + beamSpot = cms.InputTag("offlineBeamSpot"), + primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices"), + tracks = cms.InputTag("unpackedTracksAndVertices"), + minHits = cms.uint32(6), # 8 + maximumLongitudinalImpactParameter = cms.double(20), # 0.3 + minPt = cms.double(0.8), # 0.8 + maxNTracks = cms.uint32(100), # 30 + + clusterizer = cms.PSet( + seedMax3DIPSignificance = cms.double(9999.), # 9999. + seedMax3DIPValue = cms.double(9999.), # 9999. + seedMin3DIPSignificance = cms.double(1.2), # 1.2 + seedMin3DIPValue = cms.double(0.005), # 0.005 + clusterMaxDistance = cms.double(0.4), # 0.05 + clusterMaxSignificance = cms.double(4.5), # 4.5 + distanceRatio = cms.double(20), # 20 + clusterMinAngleCosine = cms.double(0.5), # 0.5 + ), + + vertexMinAngleCosine = cms.double(0.95), # 0.95 + vertexMinDLen2DSig = cms.double(2.5), # 2.5 + vertexMinDLenSig = cms.double(0.5), # 0.5 + fitterSigmacut = cms.double(3), # 3 + fitterTini = cms.double(256), # 256 + fitterRatio = cms.double(0.25), # 0.25 + useDirectVertexFitter = cms.bool(True), # True + useVertexReco = cms.bool(True), # True + vertexReco = cms.PSet( + finder = cms.string('avr'), + primcut = cms.double(1.0), # 1.0 + seccut = cms.double(3), # 3 + smoothing = cms.bool(True)) # True + ) + +displacedVertexMerger = cms.EDProducer("VertexMerger", + secondaryVertices = cms.InputTag("displacedInclusiveVertexFinder"), + maxFraction = cms.double(0.7), + minSignificance = cms.double(2)) + +# ref: https://github.com/cms-sw/cmssw/blob/2bed69b1658e4deeaef914e462741919e9183be3/RecoVertex/AdaptiveVertexFinder/plugins/VertexArbitrators.cc#L54 + +displacedTrackVertexArbitrator = cms.EDProducer("TrackVertexArbitrator", + beamSpot = cms.InputTag("offlineBeamSpot"), + primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices"), + tracks = cms.InputTag("unpackedTracksAndVertices"), + secondaryVertices = cms.InputTag("displacedVertexMerger"), + dLenFraction = cms.double(0.333), # 0.333 + dRCut = cms.double(1), # 0.4 + distCut = cms.double(0.1), # 0.04 + sigCut = cms.double(5), # 5 + fitterSigmacut = cms.double(3), # 3 + fitterTini = cms.double(256), # 256 + fitterRatio = cms.double(0.25), # 0.25 + trackMinLayers = cms.int32(4), # 4 + trackMinPt = cms.double(.4), # 0.4 + trackMinPixels = cms.int32(0) # 1 +) + +displacedInclusiveSecondaryVertices = displacedVertexMerger.clone() +displacedInclusiveSecondaryVertices.secondaryVertices = cms.InputTag("displacedTrackVertexArbitrator") +displacedInclusiveSecondaryVertices.maxFraction = 0.2 +displacedInclusiveSecondaryVertices.minSignificance = 10 + +displacedInclusiveVertexing = cms.Sequence(unpackedTracksAndVertices * displacedInclusiveVertexFinder * displacedVertexMerger * displacedTrackVertexArbitrator * displacedInclusiveSecondaryVertices) + +# MDSnano tables +from RecoMuon.MuonRechitClusterProducer.cscRechitClusterProducer_cfi import cscRechitClusterProducer +from RecoMuon.MuonRechitClusterProducer.dtRechitClusterProducer_cfi import dtRechitClusterProducer + +from PhysicsTools.NanoAOD.cscMDSshowerTable_cfi import cscMDSshowerTable +from PhysicsTools.NanoAOD.dtMDSshowerTable_cfi import dtMDSshowerTable + +cscMDSshowerTable = cscMDSshowerTable.clone( + name = cms.string("cscMDSCluster"), + recHitLabel = cms.InputTag("csc2DRecHits"), + segmentLabel = cms.InputTag("dt4DSegments"), + rpcLabel = cms.InputTag("rpcRecHits") +) + +dtMDSshowerTable = dtMDSshowerTable.clone( + name = cms.string("dtMDSCluster"), + recHitLabel = cms.InputTag("dt1DRecHits"), + rpcLabel = cms.InputTag("rpcRecHits") +) + +#DSA muon tables +from PhysicsTools.NanoAOD.simpleTrackFlatTableProducer_cfi import simpleTrackFlatTableProducer +DSAmuonsSimpleTable = simpleTrackFlatTableProducer.clone( + src="displacedStandAloneMuons", + name="DSAMuon", + doc ="Displaced standalone muon tracks", + singleton=False, + extension=False, + variables=cms.PSet( + pt = Var("pt", float, doc="pt of DSA muon"), + ptErr = Var("ptError", float, doc="ptErr of DSA muon"), + eta = Var("eta", float, doc="eta of DSA muon"), + etaErr = Var("etaError", float, doc="etaErr of DSA muon"), + phi = Var("phi", float, doc="phi of DSA muon"), + phiErr = Var("phiError", float, doc="phiErr of DSA muon"), + charge = Var("charge", float, doc="charge of DSA muon"), + dxy = Var("dxy", float, doc="dxy of DSA muon"), + dz = Var("dz", float, doc="dz of DSA muon"), + vx = Var("vx", float, doc="vx of DSA muon"), + vy = Var("vy", float, doc="vy of DSA muon"), + vz = Var("vz", float, doc="vz of DSA muon"), + chi2 = Var("chi2", float, doc="chi2 of DSA muon"), + ndof = Var("ndof", float, doc="ndof of DSA muon"), + normChi2 = Var("normalizedChi2", float, doc="normChi2 of DSA muon"), + ) +) + +DSAmuonsTable = cms.EDProducer("DSAMuonTableProducer", + name=cms.string("DSAMuon"), + dsaMuons=cms.InputTag("displacedStandAloneMuons"), + muons=cms.InputTag("linkedObjects","muons"), + primaryVertex = cms.InputTag("offlineSlimmedPrimaryVertices"), + beamspot = cms.InputTag("offlineBeamSpot") +) + +DSAmuonVertexTable = cms.EDProducer("MuonVertexTableProducer", + dsaMuons=cms.InputTag("displacedStandAloneMuons"), + patMuons=cms.InputTag("linkedObjects","muons"), + beamspot=cms.InputTag("offlineBeamSpot"), + generalTracks=cms.InputTag("generalTracks"), + primaryVertex=cms.InputTag("offlineSlimmedPrimaryVertices") +) + +from PhysicsTools.NanoAOD.simplePATMuonFlatTableProducer_cfi import simplePATMuonFlatTableProducer +PATmuonExtendedSimpleTable = simplePATMuonFlatTableProducer.clone( + src=cms.InputTag("linkedObjects","muons"), + name="Muon", + singleton=False, + extension=True, + variables=cms.PSet( + innerTrackValidFraction = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().validFraction() : -1", float, doc=""), + globalTrackNormalizedChi2 = Var("? globalTrack().isNonnull() ? globalTrack().normalizedChi2() : -1", float, doc=""), + CQChi2Position = Var("combinedQuality().chi2LocalPosition", float, doc=""), + CQTrackKink = Var("combinedQuality().trkKink", float, doc=""), + numberOfMatchedStation = Var("numberOfMatchedStations", float, doc=""), + numberOfValidPixelHits = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().hitPattern().numberOfValidPixelHits() : 0", float, doc=""), + numberOfValidTrackerHits = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().hitPattern().numberOfValidTrackerHits() : 0", float, doc=""), + trackerLayersWithMeasurement = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().hitPattern().trackerLayersWithMeasurement() : 0", float, doc=""), + numberInnerHits = Var("? globalTrack().isNonnull() ? globalTrack().hitPattern().numberOfValidMuonHits() : ? outerTrack().isNonnull() ? outerTrack().hitPattern().numberOfValidMuonHits() : 0", float, doc=""), + innerVx = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().vx() : -1", float, doc=""), + innerVy = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().vy() : -1", float, doc=""), + innerVz = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().vz() : -1", float, doc=""), + innerPt = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().pt() : -1", float, doc=""), + innerEta = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().eta() : -5", float, doc=""), + innerPhi = Var("? innerTrack().isNonnull() && innerTrack().isAvailable() ? innerTrack().phi() : -5", float, doc=""), + ) +) + +PATmuonExtendedTable = cms.EDProducer("MuonExtendedTableProducer", + name=cms.string("Muon"), + rho=cms.InputTag("fixedGridRhoFastjetAll"), + muons=cms.InputTag("linkedObjects","muons"), + primaryVertex=cms.InputTag("offlineSlimmedPrimaryVertices"), + beamspot=cms.InputTag("offlineBeamSpot"), + jets=cms.InputTag("linkedObjects","jets"), + jetsFat=cms.InputTag("slimmedJetsAK8"), + jetsSub=cms.InputTag("slimmedJetsAK8PFPuppiSoftDropPacked", "SubJets") +) + +electronVertexTable = cms.EDProducer("ElectronVertexTableProducer", + electrons=cms.InputTag("linkedObjects","electrons"), + beamspot=cms.InputTag("offlineBeamSpot"), + primaryVertex=cms.InputTag("offlineSlimmedPrimaryVertices") +) + +electronExtendedTable = cms.EDProducer("ElectronExtendedTableProducer", + name=cms.string("Electron"), + rho=cms.InputTag("fixedGridRhoFastjetAll"), + electrons=cms.InputTag("linkedObjects","electrons"), + primaryVertex=cms.InputTag("offlineSlimmedPrimaryVertices"), + jets=cms.InputTag("linkedObjects","jets"), + jetsFat=cms.InputTag("slimmedJetsAK8"), + jetsSub=cms.InputTag("slimmedJetsAK8PFPuppiSoftDropPacked", "SubJets") +) + +dispJetTable = cms.EDProducer("DispJetTableProducer", + electrons=cms.InputTag("linkedObjects","electrons"), + muons=cms.InputTag("linkedObjects","muons"), + primaryVertex = cms.InputTag("offlineSlimmedPrimaryVertices"), + secondaryVertex = cms.InputTag("displacedInclusiveSecondaryVertices") +) + +def add_dispJetTables(process): + # process.load('PhysicsTools.displacedInclusiveVertexing_cff') + process.dispJetTable = dispJetTable + process.unpackedTracksAndVertices = unpackedTracksAndVertices + process.displacedInclusiveVertexFinder = displacedInclusiveVertexFinder + process.displacedVertexMerger = displacedVertexMerger + process.displacedTrackVertexArbitrator = displacedTrackVertexArbitrator + process.displacedInclusiveSecondaryVertices = displacedInclusiveSecondaryVertices + process.dispJetTable = dispJetTable + process.dispJetTask = cms.Task(process.unpackedTracksAndVertices) + process.dispJetTask.add(process.displacedInclusiveVertexFinder) + process.dispJetTask.add(process.displacedVertexMerger) + process.dispJetTask.add(process.displacedTrackVertexArbitrator) + process.dispJetTask.add(process.displacedInclusiveSecondaryVertices) + process.dispJetTask.add(process.dispJetTable) + process.nanoTableTaskCommon.add(process.dispJetTask) + + return process + +def add_mdsTables(process): + process.ca4CSCrechitClusters = cscRechitClusterProducer + process.ca4DTrechitClusters = dtRechitClusterProducer + process.cscMDSshowerTable = cscMDSshowerTable + process.dtMDSshowerTable = dtMDSshowerTable + + process.MDSTask = cms.Task(process.ca4CSCrechitClusters) + process.MDSTask.add(process.ca4DTrechitClusters) + process.MDSTask.add(process.cscMDSshowerTable) + process.MDSTask.add(process.dtMDSshowerTable) + + process.nanoTableTaskCommon.add(process.MDSTask) + + return process + +def add_dsamuonTables(process): + process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") + process.load('TrackPropagation.SteppingHelixPropagator.SteppingHelixPropagatorAny_cfi') + + process.DSAmuonsSimpleTable = DSAmuonsSimpleTable + process.DSAmuonsTable = DSAmuonsTable + process.DSAmuonVertexTable = DSAmuonVertexTable + process.PATmuonExtendedSimpleTable = PATmuonExtendedSimpleTable + process.PATmuonExtendedTable = PATmuonExtendedTable + + process.dsamuonTask = cms.Task( + process.DSAmuonsSimpleTable, + process.DSAmuonsTable + ) + process.dsamuonVertexTask = cms.Task(process.DSAmuonVertexTable) + process.patmuonTask = cms.Task( + process.PATmuonExtendedSimpleTable, + process.PATmuonExtendedTable + ) + + process.nanoTableTaskCommon.add(process.dsamuonTask) + process.nanoTableTaskCommon.add(process.dsamuonVertexTask) + process.nanoTableTaskCommon.add(process.patmuonTask) + + return process + +def add_electronVertexTables(process): + + process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") + process.load('TrackPropagation.SteppingHelixPropagator.SteppingHelixPropagatorAny_cfi') + + process.electronVertexTable = electronVertexTable + process.electronVertexTask = cms.Task(process.electronVertexTable) + process.electronExtendedTable = electronExtendedTable + process.electronExtendedTask = cms.Task(process.electronExtendedTable) + + process.nanoTableTaskCommon.add(process.electronVertexTask) + + return process + +def add_muonExtendedTable(process): + process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") + process.load('TrackPropagation.SteppingHelixPropagator.SteppingHelixPropagatorAny_cfi') + + process.PATmuonExtendedTable = PATmuonExtendedTable + + process.muonExtendedTask = cms.Task(process.PATmuonExtendedTable) + + process.nanoTableTaskCommon.add(process.muonExtendedTask) + + return process + +def update_genParticleTable(process): + + process.genParticleTable.variables.vx = Var("vx",float, doc = "gen particle production vertex x coordinate (cm)", precision=8) + process.genParticleTable.variables.vy = Var("vy",float, doc = "gen particle production vertex y coordinate (cm)", precision=8) + process.genParticleTable.variables.vz = Var("vz",float, doc = "gen particle production vertex z coordinate (cm)", precision=8) + + process.genParticleTable.variables.px = Var("px",float, doc = "gen particle momentum x coordinate", precision=8) + process.genParticleTable.variables.py = Var("py",float, doc = "gen particle momentum y coordinate", precision=8) + process.genParticleTable.variables.pz = Var("pz",float, doc = "gen particle momentum z coordinate", precision=8) + + return process + +def add_exonanoTables(process): + + process = add_mdsTables(process) + process = add_dsamuonTables(process) + process = add_electronVertexTables(process) + process = add_dispJetTables(process) + + isMC = hasattr(process, "nanoSequenceMC") and process.schedule.contains(process.nanoSequenceMC) + + process = add_displacedtauCHSTables(process, isMC) + + if isMC: + process = update_genParticleTable(process) + + return process diff --git a/PhysicsTools/NanoAOD/test/test-exoNano.sh b/PhysicsTools/NanoAOD/test/test-exoNano.sh new file mode 100755 index 0000000000000..eba5f9bb33a33 --- /dev/null +++ b/PhysicsTools/NanoAOD/test/test-exoNano.sh @@ -0,0 +1 @@ +cmsDriver.py exoNanoMC --mc --eventcontent NANOAODSIM --datatier NANOAODSIM --conditions auto:phase1_2025_realistic --step PAT,NANO:@EXO --nThreads 4 --era Run3,Run3_2025 --filein /store/mc/RunIII2024Summer24DRPremix/TTtoLNu2Q_TuneCP5_13p6TeV_powheg-pythia8/AODSIM/140X_mcRun3_2024_realistic_v26-v2/100016/a55a360c-a746-4910-b1af-96466f009717.root -n 100 \ No newline at end of file From 91ee1bcd481b76b77848d5f4b45c530da19196f7 Mon Sep 17 00:00:00 2001 From: Lovisa Rygaard <62760237+kerstinlovisa@users.noreply.github.com> Date: Thu, 19 Feb 2026 10:48:05 +0100 Subject: [PATCH 2/7] Update PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py Co-authored-by: Matti Kortelainen --- PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py b/PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py index 7318d933da82a..79679c6839a94 100644 --- a/PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py +++ b/PhysicsTools/NanoAOD/python/custom_displacedtau_cff.py @@ -82,7 +82,7 @@ def add_displacedtauCHSTables(process, isMC): ## displaced tau part - file = "RecoTauTag/TrainingFiles/DisplacedTauId/particlenet_v1_a27159734e304ea4b7f9e0042baa9e22.pb" + file = "RecoTauTag/TrainingFiles/data/DisplacedTauId/particlenet_v1_a27159734e304ea4b7f9e0042baa9e22.pb" process.options = cms.untracked.PSet( numberOfThreads = cms.untracked.uint32(4), # Global thread count numberOfStreams = cms.untracked.uint32(4), # Should match threads From 9953a84cceaa1b663f0958a181736fc421195db9 Mon Sep 17 00:00:00 2001 From: Lovisa Date: Fri, 20 Feb 2026 10:25:06 +0100 Subject: [PATCH 3/7] code-format correction --- PhysicsTools/NanoAOD/plugins/DisTauTag.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc index f8a27692b6b96..96136acea215c 100644 --- a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc +++ b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc @@ -285,7 +285,9 @@ void DisTauTag::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& // Running inference on batch std::vector outputs; - { tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); } + { + tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); + } // Storing results for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) { From e514f2b0c2a5d7287c5ede94416a98a91871a35e Mon Sep 17 00:00:00 2001 From: Lovisa Date: Thu, 12 Mar 2026 17:46:03 +0100 Subject: [PATCH 4/7] Corrected for additional PR comments --- PhysicsTools/NanoAOD/plugins/DisTauTag.cc | 4 +- .../NanoAOD/plugins/DispJetTableProducer.cc | 91 +------------------ .../plugins/ElectronExtendedTableProducer.cc | 67 ++++++-------- .../plugins/MuonExtendedTableProducer.cc | 71 ++++++--------- .../plugins/cscMDSshowerTableProducer.cc | 4 +- PhysicsTools/NanoAOD/python/custom_exo_cff.py | 2 +- 6 files changed, 68 insertions(+), 171 deletions(-) diff --git a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc index 96136acea215c..f8a27692b6b96 100644 --- a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc +++ b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc @@ -285,9 +285,7 @@ void DisTauTag::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& // Running inference on batch std::vector outputs; - { - tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); - } + { tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); } // Storing results for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) { diff --git a/PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc b/PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc index 07972b6548991..e5f4e61a8361b 100644 --- a/PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/DispJetTableProducer.cc @@ -102,78 +102,6 @@ class DispJetTableProducer : public edm::stream::EDProducer<> { mu_IVF_tracksignedIP3Dsig; std::vector mu_IVF_signedIP2D, mu_IVF_signedIP3D, mu_IVF_signedIP2Dsig, mu_IVF_signedIP3Dsig; - el_idx.clear(); - el_lIVF_match.clear(); - el_IVF_df.clear(); - el_IVF_ntracks.clear(); - el_IVF_elid.clear(); - el_IVF_x.clear(); - el_IVF_y.clear(); - el_IVF_z.clear(); - el_IVF_cx.clear(); - el_IVF_cy.clear(); - el_IVF_cz.clear(); - el_IVF_chi2.clear(); - el_IVF_pt.clear(); - el_IVF_eta.clear(); - el_IVF_phi.clear(); - el_IVF_E.clear(); - el_IVF_mass.clear(); - el_IVF_trackcharge.clear(); - el_IVF_trackelid.clear(); - el_IVF_trackvtxid.clear(); - el_IVF_trackpt.clear(); - el_IVF_tracketa.clear(); - el_IVF_trackphi.clear(); - el_IVF_trackE.clear(); - el_IVF_trackdxy.clear(); - el_IVF_trackdz.clear(); - el_IVF_tracksignedIP2D.clear(); - el_IVF_tracksignedIP2Dsig.clear(); - el_IVF_tracksignedIP3D.clear(); - el_IVF_tracksignedIP3Dsig.clear(); - - el_IVF_signedIP2D.clear(); - el_IVF_signedIP2Dsig.clear(); - el_IVF_signedIP3D.clear(); - el_IVF_signedIP3Dsig.clear(); - - mu_idx.clear(); - mu_lIVF_match.clear(); - mu_IVF_df.clear(); - mu_IVF_ntracks.clear(); - mu_IVF_muid.clear(); - mu_IVF_x.clear(); - mu_IVF_y.clear(); - mu_IVF_z.clear(); - mu_IVF_cx.clear(); - mu_IVF_cy.clear(); - mu_IVF_cz.clear(); - mu_IVF_chi2.clear(); - mu_IVF_pt.clear(); - mu_IVF_eta.clear(); - mu_IVF_phi.clear(); - mu_IVF_E.clear(); - mu_IVF_mass.clear(); - mu_IVF_trackcharge.clear(); - mu_IVF_trackmuid.clear(); - mu_IVF_trackvtxid.clear(); - mu_IVF_trackpt.clear(); - mu_IVF_tracketa.clear(); - mu_IVF_trackphi.clear(); - mu_IVF_trackE.clear(); - mu_IVF_trackdxy.clear(); - mu_IVF_trackdz.clear(); - mu_IVF_tracksignedIP2D.clear(); - mu_IVF_tracksignedIP2Dsig.clear(); - mu_IVF_tracksignedIP3D.clear(); - mu_IVF_tracksignedIP3Dsig.clear(); - - mu_IVF_signedIP2D.clear(); - mu_IVF_signedIP2Dsig.clear(); - mu_IVF_signedIP3D.clear(); - mu_IVF_signedIP3Dsig.clear(); - int ntrack_max = 100; int nElectronsSel = 0; int nMuonsSel = 0; @@ -181,12 +109,9 @@ class DispJetTableProducer : public edm::stream::EDProducer<> { for (unsigned int i = 0; i < nElectrons; i++) { const pat::Electron& el = electrons[i]; - if (el.gsfTrack().isNull()) - continue; - if (el.pt() < 7) - continue; - if (fabs(el.eta()) > 2.5) + if (el.gsfTrack().isNull() || el.pt() < 7 || fabs(el.eta()) > 2.5) { continue; + } el_idx.push_back(i); el_lIVF_match.push_back(false); @@ -284,16 +209,10 @@ class DispJetTableProducer : public edm::stream::EDProducer<> { for (unsigned int i = 0; i < nMuons; i++) { const pat::Muon& mu = muons[i]; - if (mu.innerTrack().isNull()) - continue; - if (mu.pt() < 5) - continue; - if (fabs(mu.eta()) > 2.4) - continue; - if (!mu.isPFMuon()) - continue; - if (!(mu.isTrackerMuon() || mu.isGlobalMuon())) + if (mu.innerTrack().isNull() || mu.pt() < 5 || fabs(mu.eta()) > 2.4 || !mu.isPFMuon() || + !(mu.isTrackerMuon() || mu.isGlobalMuon())) { continue; + } mu_idx.push_back(i); mu_lIVF_match.push_back(false); diff --git a/PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc b/PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc index 61ee491eb4a77..a9bb60f4eedd8 100644 --- a/PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/ElectronExtendedTableProducer.cc @@ -60,11 +60,11 @@ class ElectronExtendedTableProducer : public edm::global::EDProducer<> { const std::vector& jets, const reco::Vertex& vertex, const double rho, - std::vector* jetIdx, + std::vector& jetIdx, std::vector relIso0p4, - std::vector* jetPtRatio, - std::vector* jetPtRel, - std::vector* jrtSelectedChargedMultiplicity) const; + std::vector& jetPtRatio, + std::vector& jetPtRel, + std::vector& jrtSelectedChargedMultiplicity) const; std::string name_; edm::EDGetTokenT rhoTag_; @@ -135,7 +135,7 @@ void ElectronExtendedTableProducer::produce(edm::StreamID, edm::Event& iEvent, c relIso0p4.push_back(getPFIso(electron)); fillLeptonJetVariables( - &electron, jets, pv, rho, &jetIdx, relIso0p4, &jetPtRatio, &jetPtRel, &jetSelectedChargedMultiplicity); + &electron, jets, pv, rho, jetIdx, relIso0p4, jetPtRatio, jetPtRel, jetSelectedChargedMultiplicity); const reco::Candidate* el_cand = dynamic_cast(&electron); jetFatIdx.push_back(findMatchedJet(*el_cand, jetFat)); @@ -217,21 +217,21 @@ void ElectronExtendedTableProducer::fillLeptonJetVariables(const reco::GsfElectr const std::vector& jets, const reco::Vertex& vertex, const double rho, - std::vector* jetIdx, + std::vector& jetIdx, std::vector relIso0p4, - std::vector* jetPtRatio, - std::vector* jetPtRel, - std::vector* jetSelectedChargedMultiplicity) const { + std::vector& jetPtRatio, + std::vector& jetPtRel, + std::vector& jetSelectedChargedMultiplicity) const { const reco::Candidate* cand = dynamic_cast(el); int matchedJetIdx = findMatchedJet(*cand, jets); - jetIdx->push_back(matchedJetIdx); + jetIdx.push_back(matchedJetIdx); if (matchedJetIdx < 0) { float ptRatio = (1. / (1. + relIso0p4.back())); - jetPtRatio->push_back(ptRatio); - jetPtRel->push_back(0); - jetSelectedChargedMultiplicity->push_back(0); + jetPtRatio.push_back(ptRatio); + jetPtRel.push_back(0); + jetSelectedChargedMultiplicity.push_back(0); } else { const pat::Jet& jet = jets[matchedJetIdx]; auto rawJetP4 = jet.correctedP4("Uncorrected"); @@ -240,9 +240,9 @@ void ElectronExtendedTableProducer::fillLeptonJetVariables(const reco::GsfElectr bool leptonEqualsJet = ((rawJetP4 - leptonP4).P() < 1e-4); if (leptonEqualsJet) { - jetPtRatio->push_back(1); - jetPtRel->push_back(0); - jetSelectedChargedMultiplicity->push_back(0); + jetPtRatio.push_back(1); + jetPtRel.push_back(0); + jetSelectedChargedMultiplicity.push_back(0); } else { auto L1JetP4 = jet.correctedP4("L1FastJet"); double L2L3JEC = jet.pt() / L1JetP4.pt(); @@ -250,37 +250,28 @@ void ElectronExtendedTableProducer::fillLeptonJetVariables(const reco::GsfElectr float ptRatio = cand->pt() / lepAwareJetP4.pt(); float ptRel = leptonP4.Vect().Cross((lepAwareJetP4 - leptonP4).Vect().Unit()).R(); - jetPtRatio->push_back(ptRatio); - jetPtRel->push_back(ptRel); - jetSelectedChargedMultiplicity->push_back(0); + jetPtRatio.push_back(ptRatio); + jetPtRel.push_back(ptRel); + jetSelectedChargedMultiplicity.push_back(0); for (const auto& daughterPtr : jet.daughterPtrVector()) { const pat::PackedCandidate& daughter = *((const pat::PackedCandidate*)daughterPtr.get()); - if (daughter.charge() == 0) - continue; - if (daughter.fromPV() < 2) - continue; - if (reco::deltaR(daughter, *cand) > 0.4) - continue; - if (!daughter.hasTrackDetails()) + if (daughter.charge() == 0 || daughter.fromPV() < 2 || reco::deltaR(daughter, *cand) > 0.4 || + !daughter.hasTrackDetails()) { continue; + } auto daughterTrack = daughter.pseudoTrack(); - if (daughterTrack.pt() <= 1) - continue; - if (daughterTrack.hitPattern().numberOfValidHits() < 8) - continue; - if (daughterTrack.hitPattern().numberOfValidPixelHits() < 2) + if (daughterTrack.pt() <= 1 || daughterTrack.hitPattern().numberOfValidHits() < 8 || + daughterTrack.hitPattern().numberOfValidPixelHits() < 2 || daughterTrack.normalizedChi2() >= 5 || + std::abs(daughterTrack.dz(vertex.position())) >= 17 || + std::abs(daughterTrack.dxy(vertex.position())) >= 0.2) { continue; - if (daughterTrack.normalizedChi2() >= 5) - continue; - if (std::abs(daughterTrack.dz(vertex.position())) >= 17) - continue; - if (std::abs(daughterTrack.dxy(vertex.position())) >= 0.2) - continue; - ++jetSelectedChargedMultiplicity->back(); + } + + ++jetSelectedChargedMultiplicity.back(); } } } diff --git a/PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc b/PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc index 18d12db42cb74..fb4653eac2ba8 100644 --- a/PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/MuonExtendedTableProducer.cc @@ -63,11 +63,11 @@ class MuonExtendedTableProducer : public edm::global::EDProducer<> { const std::vector& jets, const reco::Vertex& vertex, const double rho, - std::vector* jetIdx, + std::vector& jetIdx, std::vector relIso0p4, - std::vector* jetPtRatio, - std::vector* jetPtRel, - std::vector* jetSelectedChargedMultiplicity) const; + std::vector& jetPtRatio, + std::vector& jetPtRel, + std::vector& jetSelectedChargedMultiplicity) const; std::string name_; edm::EDGetTokenT rhoTag_; @@ -144,7 +144,7 @@ void MuonExtendedTableProducer::produce(edm::StreamID, edm::Event& iEvent, const relIso0p4.push_back(getPFIso(muon)); fillLeptonJetVariables( - &muon, jets, pv, rho, &jetIdx, relIso0p4, &jetPtRatio, &jetPtRel, &jetSelectedChargedMultiplicity); + &muon, jets, pv, rho, jetIdx, relIso0p4, jetPtRatio, jetPtRel, jetSelectedChargedMultiplicity); const reco::Candidate* mu_cand = dynamic_cast(&muon); jetFatIdx.push_back(findMatchedJet(*mu_cand, jetFat)); @@ -246,8 +246,6 @@ bool isSourceCandidatePtrMatch(const T1& lhs, const T2& rhs) { } int MuonExtendedTableProducer::findMatchedJet(const reco::Candidate& lepton, const std::vector& jets) const { - int iJet = -1; - unsigned int nJets = jets.size(); for (unsigned int i = 0; i < nJets; i++) { @@ -257,28 +255,28 @@ int MuonExtendedTableProducer::findMatchedJet(const reco::Candidate& lepton, con } } - return iJet; + return -1; } void MuonExtendedTableProducer::fillLeptonJetVariables(const reco::Muon* mu, const std::vector& jets, const reco::Vertex& vertex, const double rho, - std::vector* jetIdx, + std::vector& jetIdx, std::vector relIso0p4, - std::vector* jetPtRatio, - std::vector* jetPtRel, - std::vector* jetSelectedChargedMultiplicity) const { + std::vector& jetPtRatio, + std::vector& jetPtRel, + std::vector& jetSelectedChargedMultiplicity) const { const reco::Candidate* cand = dynamic_cast(mu); int matchedJetIdx = findMatchedJet(*cand, jets); - jetIdx->push_back(matchedJetIdx); + jetIdx.push_back(matchedJetIdx); if (matchedJetIdx < 0) { float ptRatio = (1. / (1. + relIso0p4.back())); - jetPtRatio->push_back(ptRatio); - jetPtRel->push_back(0); - jetSelectedChargedMultiplicity->push_back(0); + jetPtRatio.push_back(ptRatio); + jetPtRel.push_back(0); + jetSelectedChargedMultiplicity.push_back(0); } else { const pat::Jet& jet = jets[matchedJetIdx]; auto rawJetP4 = jet.correctedP4("Uncorrected"); @@ -287,9 +285,9 @@ void MuonExtendedTableProducer::fillLeptonJetVariables(const reco::Muon* mu, bool leptonEqualsJet = ((rawJetP4 - leptonP4).P() < 1e-4); if (leptonEqualsJet) { - jetPtRatio->push_back(1); - jetPtRel->push_back(0); - jetSelectedChargedMultiplicity->push_back(0); + jetPtRatio.push_back(1); + jetPtRel.push_back(0); + jetSelectedChargedMultiplicity.push_back(0); } else { auto L1JetP4 = jet.correctedP4("L1FastJet"); double L2L3JEC = jet.pt() / L1JetP4.pt(); @@ -297,37 +295,28 @@ void MuonExtendedTableProducer::fillLeptonJetVariables(const reco::Muon* mu, float ptRatio = cand->pt() / lepAwareJetP4.pt(); float ptRel = leptonP4.Vect().Cross((lepAwareJetP4 - leptonP4).Vect().Unit()).R(); - jetPtRatio->push_back(ptRatio); - jetPtRel->push_back(ptRel); - jetSelectedChargedMultiplicity->push_back(0); + jetPtRatio.push_back(ptRatio); + jetPtRel.push_back(ptRel); + jetSelectedChargedMultiplicity.push_back(0); for (const auto& daughterPtr : jet.daughterPtrVector()) { const pat::PackedCandidate& daughter = *((const pat::PackedCandidate*)daughterPtr.get()); - if (daughter.charge() == 0) - continue; - if (daughter.fromPV() < 2) - continue; - if (reco::deltaR(daughter, *cand) > 0.4) - continue; - if (!daughter.hasTrackDetails()) + if (daughter.charge() == 0 || daughter.fromPV() < 2 || reco::deltaR(daughter, *cand) > 0.4 || + !daughter.hasTrackDetails()) { continue; + } auto daughterTrack = daughter.pseudoTrack(); - if (daughterTrack.pt() <= 1) - continue; - if (daughterTrack.hitPattern().numberOfValidHits() < 8) + if (daughterTrack.pt() <= 1 || daughterTrack.hitPattern().numberOfValidHits() < 8 || + daughterTrack.hitPattern().numberOfValidPixelHits() < 2 || daughterTrack.normalizedChi2() >= 5 || + std::abs(daughterTrack.dz(vertex.position())) >= 17 || + std::abs(daughterTrack.dxy(vertex.position())) >= 0.2) { continue; - if (daughterTrack.hitPattern().numberOfValidPixelHits() < 2) - continue; - if (daughterTrack.normalizedChi2() >= 5) - continue; - if (std::abs(daughterTrack.dz(vertex.position())) >= 17) - continue; - if (std::abs(daughterTrack.dxy(vertex.position())) >= 0.2) - continue; - ++jetSelectedChargedMultiplicity->back(); + } + + ++jetSelectedChargedMultiplicity.back(); } } } diff --git a/PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc b/PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc index 014b5fb433d0f..b1d8316cb61a6 100644 --- a/PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc +++ b/PhysicsTools/NanoAOD/plugins/cscMDSshowerTableProducer.cc @@ -40,7 +40,7 @@ class cscMDSshowerTableProducer : public edm::global::EDProducer<> { dtgeometryToken_(esConsumes()), rpcgeometryToken_(esConsumes()), inputToken_(consumes(iConfig.getParameter("recHitLabel"))), - dtSegmentToken_(consumes(iConfig.getParameter("segmentLabel"))), + dtSegmentToken_(consumes(iConfig.getParameter("dtSegmentLabel"))), rpchitToken_(consumes(iConfig.getParameter("rpcLabel"))), rParam_(iConfig.getParameter("rParam")), nRechitMin_(iConfig.getParameter("nRechitMin")), @@ -58,7 +58,7 @@ class cscMDSshowerTableProducer : public edm::global::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add("recHitLabel")->setComment("input cscRechit collection"); - desc.add("segmentLabel")->setComment("input dt segment collection for veto"); + desc.add("dtSegmentLabel")->setComment("input dt segment collection for veto"); desc.add("rpcLabel")->setComment("input rpcRechit collection for veto"); desc.add("rParam", 0.4); desc.add("nRechitMin", 50); diff --git a/PhysicsTools/NanoAOD/python/custom_exo_cff.py b/PhysicsTools/NanoAOD/python/custom_exo_cff.py index 9cedea34d92a1..f356653c52035 100644 --- a/PhysicsTools/NanoAOD/python/custom_exo_cff.py +++ b/PhysicsTools/NanoAOD/python/custom_exo_cff.py @@ -88,7 +88,7 @@ cscMDSshowerTable = cscMDSshowerTable.clone( name = cms.string("cscMDSCluster"), recHitLabel = cms.InputTag("csc2DRecHits"), - segmentLabel = cms.InputTag("dt4DSegments"), + dtSegmentLabel = cms.InputTag("dt4DSegments"), rpcLabel = cms.InputTag("rpcRecHits") ) From 8263cc138f273cb8162f085e21439f934f8309b0 Mon Sep 17 00:00:00 2001 From: Lovisa Date: Thu, 12 Mar 2026 17:53:41 +0100 Subject: [PATCH 5/7] code format --- PhysicsTools/NanoAOD/plugins/DisTauTag.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc index f8a27692b6b96..96136acea215c 100644 --- a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc +++ b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc @@ -285,7 +285,9 @@ void DisTauTag::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& // Running inference on batch std::vector outputs; - { tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); } + { + tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); + } // Storing results for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) { From abeed7be932b6318f11d46f42774bf3493cbb063 Mon Sep 17 00:00:00 2001 From: Lovisa Date: Tue, 7 Apr 2026 08:37:05 +0200 Subject: [PATCH 6/7] Updating DisTauTag.cc test_vector to test_score --- PhysicsTools/NanoAOD/plugins/DisTauTag.cc | 25 ++++++++--------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc index 96136acea215c..5bb8099863e42 100644 --- a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc +++ b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc @@ -32,16 +32,12 @@ #include "PhysicsTools/NanoAOD/interface/DisTauTagScaling.h" -void test_vector(std::vector& values) { - for (auto& value : values) { - if (std::isnan(value)) { - throw std::runtime_error("DisTauTag score output: NaN detected."); - } else if (std::isinf(value)) { - throw std::runtime_error("DisTauTag score output: Infinity detected."); - } else if (!std::isfinite(value)) { - throw std::runtime_error("DisTauTag score output: Non-standard value detected."); - } +std::optional test_score(float value) { + if (!std::isfinite(value)) { + edm::LogWarning("DisTauTag") << "Non-finite score detected (" << value << "); dropping value."; + return std::nullopt; } + return value; } class DisTauTag : public edm::global::EDProducer<> { @@ -285,15 +281,13 @@ void DisTauTag::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& // Running inference on batch std::vector outputs; - { - tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); - } + { tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); } // Storing results for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) { const size_t jetIndex = jet_start + batch_jet_idx; - v_score0[jetIndex] = outputs[0].matrix()(batch_jet_idx, 0); - v_score1[jetIndex] = outputs[0].matrix()(batch_jet_idx, 1); + v_score0[jetIndex] = test_score(outputs[0].matrix()(batch_jet_idx, 0)).value_or(-9.f); + v_score1[jetIndex] = test_score(outputs[0].matrix()(batch_jet_idx, 1)).value_or(-9.f); } // Save inputs if requested (for debugging) @@ -321,9 +315,6 @@ void DisTauTag::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& } } - test_vector(v_score0); - test_vector(v_score1); - std::unique_ptr> vm_score0(new edm::ValueMap()); edm::ValueMap::Filler filler_score0(*vm_score0); filler_score0.insert(jets, v_score0.begin(), v_score0.end()); From e890a13540bda295bf0d6233e763caa02d875285 Mon Sep 17 00:00:00 2001 From: Lovisa Date: Tue, 7 Apr 2026 08:43:55 +0200 Subject: [PATCH 7/7] code format fix --- PhysicsTools/NanoAOD/plugins/DisTauTag.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc index 5bb8099863e42..92d5f0d0ea3c5 100644 --- a/PhysicsTools/NanoAOD/plugins/DisTauTag.cc +++ b/PhysicsTools/NanoAOD/plugins/DisTauTag.cc @@ -281,7 +281,9 @@ void DisTauTag::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& // Running inference on batch std::vector outputs; - { tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); } + { + tensorflow::run(session_, {{"input_1", input_1}, {"input_2", input_2}}, {"final_out"}, &outputs); + } // Storing results for (size_t batch_jet_idx = 0; batch_jet_idx < current_batch_size; ++batch_jet_idx) {