File indexing completed on 2024-10-08 23:09:56
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021
0022
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/stream/EDProducer.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include "RecoJets/JetPlusTracks/plugins/JetPlusTrackProducer.h"
0031 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0032 #include "DataFormats/JetReco/interface/CaloJet.h"
0033 #include "DataFormats/JetReco/interface/JPTJetCollection.h"
0034 #include "DataFormats/JetReco/interface/JPTJet.h"
0035 #include "DataFormats/JetReco/interface/TrackJetCollection.h"
0036 #include "DataFormats/JetReco/interface/TrackJet.h"
0037 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0038 #include "DataFormats/TrackReco/interface/Track.h"
0039 #include "DataFormats/JetReco/interface/Jet.h"
0040 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0041 #include "DataFormats/VertexReco/interface/Vertex.h"
0042 #include "DataFormats/Math/interface/deltaPhi.h"
0043 #include "DataFormats/Math/interface/deltaR.h"
0044 #include <string>
0045
0046 using namespace std;
0047 using namespace jpt;
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060 JetPlusTrackProducer::JetPlusTrackProducer(const edm::ParameterSet& iConfig) {
0061
0062 src_ = iConfig.getParameter<edm::InputTag>("src");
0063 srcTrackJets_ = iConfig.getParameter<edm::InputTag>("srcTrackJets");
0064 alias_ = iConfig.getUntrackedParameter<string>("alias");
0065 srcPVs_ = iConfig.getParameter<edm::InputTag>("srcPVs");
0066 vectorial_ = iConfig.getParameter<bool>("VectorialCorrection");
0067 useZSP_ = iConfig.getParameter<bool>("UseZSP");
0068 ptCUT_ = iConfig.getParameter<double>("ptCUT");
0069 dRcone_ = iConfig.getParameter<double>("dRcone");
0070 usePAT_ = iConfig.getParameter<bool>("UsePAT");
0071
0072 mJPTalgo = new JetPlusTrackCorrector(iConfig, consumesCollector());
0073 if (useZSP_)
0074 mZSPalgo = new ZSPJPTJetCorrector(iConfig);
0075
0076 produces<reco::JPTJetCollection>().setBranchAlias(alias_);
0077 produces<reco::CaloJetCollection>().setBranchAlias("ak4CaloJetsJPT");
0078
0079 input_jets_token_ = consumes<edm::View<reco::CaloJet> >(src_);
0080 input_addjets_token_ = consumes<edm::View<reco::CaloJet> >(iConfig.getParameter<edm::InputTag>("srcAddCaloJets"));
0081 input_trackjets_token_ = consumes<edm::View<reco::TrackJet> >(srcTrackJets_);
0082 input_vertex_token_ = consumes<reco::VertexCollection>(srcPVs_);
0083 mExtrapolations_ =
0084 consumes<std::vector<reco::TrackExtrapolation> >(iConfig.getParameter<edm::InputTag>("extrapolations"));
0085 }
0086
0087 JetPlusTrackProducer::~JetPlusTrackProducer() {
0088
0089
0090 }
0091
0092
0093
0094
0095 bool sort_by_pt(const reco::JPTJet& a, const reco::JPTJet& b) { return (a.pt() > b.pt()); }
0096
0097
0098 void JetPlusTrackProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0099 using namespace edm;
0100
0101 auto const& jets_h = iEvent.get(input_jets_token_);
0102 auto const& addjets_h = iEvent.get(input_addjets_token_);
0103 auto const& iExtrapolations = iEvent.get(mExtrapolations_);
0104 edm::RefProd<reco::CaloJetCollection> pOut1RefProd = iEvent.getRefBeforePut<reco::CaloJetCollection>();
0105 edm::Ref<reco::CaloJetCollection>::key_type idxCaloJet = 0;
0106
0107 auto pOut = std::make_unique<reco::JPTJetCollection>();
0108 auto pOut1 = std::make_unique<reco::CaloJetCollection>();
0109
0110 for (auto const& jet : iEvent.get(input_trackjets_token_)) {
0111 int icalo = -1;
0112 int i = 0;
0113 for (auto const& oldjet : addjets_h) {
0114 double dr2 = deltaR2(jet, oldjet);
0115 if (dr2 <= dRcone_ * dRcone_) {
0116 icalo = i;
0117 }
0118 i++;
0119 }
0120 if (icalo < 0)
0121 continue;
0122 auto const& mycalo = addjets_h[icalo];
0123 std::vector<edm::Ptr<reco::Track> > tracksinjet = jet.tracks();
0124 reco::TrackRefVector tracksincalo;
0125 reco::TrackRefVector tracksinvert;
0126 for (auto const& itrack : tracksinjet) {
0127 for (auto const& ixtrp : iExtrapolations) {
0128 if (ixtrp.positions().empty())
0129 continue;
0130 if (usePAT_) {
0131 double mydphi = deltaPhi(ixtrp.track()->phi(), itrack->phi());
0132 if (fabs(ixtrp.track()->pt() - itrack->pt()) > 0.001 || fabs(ixtrp.track()->eta() - itrack->eta()) > 0.001 ||
0133 mydphi > 0.001)
0134 continue;
0135 } else {
0136 if (itrack.id() != ixtrp.track().id() || itrack.key() != ixtrp.track().key())
0137 continue;
0138 }
0139 tracksinvert.push_back(ixtrp.track());
0140 reco::TrackBase::Point const& point = ixtrp.positions().at(0);
0141 double dr2 = deltaR2(jet, point);
0142 if (dr2 <= dRcone_ * dRcone_) {
0143 tracksincalo.push_back(ixtrp.track());
0144 }
0145 }
0146 }
0147
0148 const reco::TrackJet& corrected = jet;
0149 math::XYZTLorentzVector p4;
0150 jpt::MatchedTracks pions;
0151 jpt::MatchedTracks muons;
0152 jpt::MatchedTracks elecs;
0153
0154 mJPTalgo->correction(corrected, mycalo, iEvent, iSetup, tracksinvert, tracksincalo, p4, pions, muons, elecs);
0155 if (p4.pt() > ptCUT_) {
0156 reco::JPTJet::Specific jptspe;
0157 jptspe.pionsInVertexInCalo = pions.inVertexInCalo_;
0158 jptspe.pionsInVertexOutCalo = pions.inVertexOutOfCalo_;
0159 jptspe.pionsOutVertexInCalo = pions.outOfVertexInCalo_;
0160 jptspe.muonsInVertexInCalo = muons.inVertexInCalo_;
0161 jptspe.muonsInVertexOutCalo = muons.inVertexOutOfCalo_;
0162 jptspe.muonsOutVertexInCalo = muons.outOfVertexInCalo_;
0163 jptspe.elecsInVertexInCalo = elecs.inVertexInCalo_;
0164 jptspe.elecsInVertexOutCalo = elecs.inVertexOutOfCalo_;
0165 jptspe.elecsOutVertexInCalo = elecs.outOfVertexInCalo_;
0166 reco::CaloJetRef myjet(pOut1RefProd, idxCaloJet++);
0167 jptspe.theCaloJetRef = edm::RefToBase<reco::Jet>(myjet);
0168 jptspe.JPTSeed = 1;
0169 reco::JPTJet fJet(p4, jet.primaryVertex()->position(), jptspe, mycalo.getJetConstituents());
0170 pOut->push_back(fJet);
0171 pOut1->push_back(mycalo);
0172 }
0173 }
0174
0175 int iJet = 0;
0176 for (auto const& oldjet : jets_h) {
0177 reco::CaloJet corrected = oldjet;
0178
0179
0180 double factorZSP = 1.;
0181 if (useZSP_)
0182 factorZSP = mZSPalgo->correction(corrected, iEvent, iSetup);
0183 corrected.scaleEnergy(factorZSP);
0184
0185 math::XYZTLorentzVector p4;
0186
0187 jpt::MatchedTracks pions;
0188 jpt::MatchedTracks muons;
0189 jpt::MatchedTracks elecs;
0190 bool validMatches = false;
0191
0192 if (!vectorial_) {
0193 double scaleJPT = mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, pions, muons, elecs, validMatches);
0194 p4 = math::XYZTLorentzVector(corrected.px() * scaleJPT,
0195 corrected.py() * scaleJPT,
0196 corrected.pz() * scaleJPT,
0197 corrected.energy() * scaleJPT);
0198 } else {
0199 mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, p4, pions, muons, elecs, validMatches);
0200 }
0201
0202 reco::JPTJet::Specific specific;
0203
0204 if (validMatches) {
0205 specific.pionsInVertexInCalo = pions.inVertexInCalo_;
0206 specific.pionsInVertexOutCalo = pions.inVertexOutOfCalo_;
0207 specific.pionsOutVertexInCalo = pions.outOfVertexInCalo_;
0208 specific.muonsInVertexInCalo = muons.inVertexInCalo_;
0209 specific.muonsInVertexOutCalo = muons.inVertexOutOfCalo_;
0210 specific.muonsOutVertexInCalo = muons.outOfVertexInCalo_;
0211 specific.elecsInVertexInCalo = elecs.inVertexInCalo_;
0212 specific.elecsInVertexOutCalo = elecs.inVertexOutOfCalo_;
0213 specific.elecsOutVertexInCalo = elecs.outOfVertexInCalo_;
0214 }
0215
0216
0217 specific.theCaloJetRef = edm::RefToBase<reco::Jet>(jets_h.refAt(iJet));
0218 specific.mResponseOfChargedWithEff = (float)mJPTalgo->getResponseOfChargedWithEff();
0219 specific.mResponseOfChargedWithoutEff = (float)mJPTalgo->getResponseOfChargedWithoutEff();
0220 specific.mSumPtOfChargedWithEff = (float)mJPTalgo->getSumPtWithEff();
0221 specific.mSumPtOfChargedWithoutEff = (float)mJPTalgo->getSumPtWithoutEff();
0222 specific.mSumEnergyOfChargedWithEff = (float)mJPTalgo->getSumEnergyWithEff();
0223 specific.mSumEnergyOfChargedWithoutEff = (float)mJPTalgo->getSumEnergyWithoutEff();
0224 specific.mChargedHadronEnergy = (float)mJPTalgo->getSumEnergyWithoutEff();
0225
0226
0227
0228 double deR2Tr = 0.;
0229 double deEta2Tr = 0.;
0230 double dePhi2Tr = 0.;
0231 double Zch = 0.;
0232 double Pout2 = 0.;
0233 double Pout = 0.;
0234 double denominator_tracks = 0.;
0235 int ntracks = 0;
0236
0237 for (reco::TrackRefVector::const_iterator it = pions.inVertexInCalo_.begin(); it != pions.inVertexInCalo_.end();
0238 it++) {
0239 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0240 double deEta = (*it)->eta() - p4.eta();
0241 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0242 if ((**it).ptError() / (**it).pt() < 0.1) {
0243 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0244 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0245 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0246 denominator_tracks = denominator_tracks + (*it)->pt();
0247 Zch = Zch + (*it)->pt();
0248
0249 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0250 ntracks++;
0251 }
0252 }
0253 for (reco::TrackRefVector::const_iterator it = muons.inVertexInCalo_.begin(); it != muons.inVertexInCalo_.end();
0254 it++) {
0255 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0256 double deEta = (*it)->eta() - p4.eta();
0257 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0258 if ((**it).ptError() / (**it).pt() < 0.1) {
0259 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0260 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0261 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0262 denominator_tracks = denominator_tracks + (*it)->pt();
0263 Zch = Zch + (*it)->pt();
0264
0265 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0266 ntracks++;
0267 }
0268 }
0269 for (reco::TrackRefVector::const_iterator it = elecs.inVertexInCalo_.begin(); it != elecs.inVertexInCalo_.end();
0270 it++) {
0271 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0272 double deEta = (*it)->eta() - p4.eta();
0273 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0274 if ((**it).ptError() / (**it).pt() < 0.1) {
0275 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0276 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0277 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0278 denominator_tracks = denominator_tracks + (*it)->pt();
0279 Zch = Zch + (*it)->pt();
0280
0281 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0282 ntracks++;
0283 }
0284 }
0285
0286 for (reco::TrackRefVector::const_iterator it = pions.inVertexOutOfCalo_.begin();
0287 it != pions.inVertexOutOfCalo_.end();
0288 it++) {
0289 Zch = Zch + (*it)->pt();
0290 }
0291 for (reco::TrackRefVector::const_iterator it = muons.inVertexOutOfCalo_.begin();
0292 it != muons.inVertexOutOfCalo_.end();
0293 it++) {
0294 Zch = Zch + (*it)->pt();
0295 }
0296 for (reco::TrackRefVector::const_iterator it = elecs.inVertexOutOfCalo_.begin();
0297 it != elecs.inVertexOutOfCalo_.end();
0298 it++) {
0299 Zch = Zch + (*it)->pt();
0300 }
0301
0302 if (mJPTalgo->getSumPtForBeta() > 0.)
0303 Zch = Zch / mJPTalgo->getSumPtForBeta();
0304
0305 if (ntracks > 0) {
0306 Pout = sqrt(fabs(Pout2)) / ntracks;
0307 }
0308 if (denominator_tracks != 0) {
0309 deR2Tr = deR2Tr / denominator_tracks;
0310 deEta2Tr = deEta2Tr / denominator_tracks;
0311 dePhi2Tr = dePhi2Tr / denominator_tracks;
0312 }
0313
0314 specific.R2momtr = deR2Tr;
0315 specific.Eta2momtr = deEta2Tr;
0316 specific.Phi2momtr = dePhi2Tr;
0317 specific.Pout = Pout;
0318 specific.Zch = Zch;
0319
0320
0321 reco::Particle::Point vertex_ = reco::Jet::Point(0, 0, 0);
0322
0323
0324 edm::Handle<reco::VertexCollection> pvCollection;
0325 iEvent.getByToken(input_vertex_token_, pvCollection);
0326 if (pvCollection.isValid() && !pvCollection->empty())
0327 vertex_ = pvCollection->begin()->position();
0328
0329 reco::JPTJet fJet(p4, vertex_, specific, corrected.getJetConstituents());
0330 iJet++;
0331
0332
0333 if (fJet.pt() > ptCUT_)
0334 pOut->push_back(fJet);
0335 }
0336 std::sort(pOut->begin(), pOut->end(), sort_by_pt);
0337 iEvent.put(std::move(pOut1));
0338 iEvent.put(std::move(pOut));
0339 }
0340
0341
0342