File indexing completed on 2023-03-17 11:18:29
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 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029
0030 #include "RecoJets/JetPlusTracks/plugins/JetPlusTrackProducerAA.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/TrackReco/interface/TrackFwd.h"
0036 #include "DataFormats/TrackReco/interface/Track.h"
0037 #include "DataFormats/JetReco/interface/Jet.h"
0038 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0039 #include "DataFormats/VertexReco/interface/Vertex.h"
0040 #include "DataFormats/Math/interface/deltaPhi.h"
0041 #include "DataFormats/Math/interface/deltaR.h"
0042
0043
0044 #include "RecoJets/JetAssociationAlgorithms/interface/JetTracksAssociationXtrpCalo.h"
0045 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0046 #include "DataFormats/GeometrySurface/interface/Plane.h"
0047 #include "DataFormats/Math/interface/deltaR.h"
0048 #include "DataFormats/Math/interface/Vector3D.h"
0049 #include "MagneticField/Engine/interface/MagneticField.h"
0050 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0051 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0052 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0053 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0054 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0055 #include "DataFormats/DetId/interface/DetId.h"
0056 #include "DataFormats/JetReco/interface/CaloJet.h"
0057 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0058
0059
0060 #include <string>
0061
0062 using namespace std;
0063 using namespace jpt;
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076 JetPlusTrackProducerAA::JetPlusTrackProducerAA(const edm::ParameterSet& iConfig) {
0077
0078 src = iConfig.getParameter<edm::InputTag>("src");
0079 alias = iConfig.getUntrackedParameter<string>("alias");
0080 mTracks = iConfig.getParameter<edm::InputTag>("tracks");
0081 srcPVs_ = iConfig.getParameter<edm::InputTag>("srcPVs");
0082 vectorial_ = iConfig.getParameter<bool>("VectorialCorrection");
0083 useZSP = iConfig.getParameter<bool>("UseZSP");
0084 std::string tq = iConfig.getParameter<std::string>("TrackQuality");
0085 trackQuality_ = reco::TrackBase::qualityByName(tq);
0086 mConeSize = iConfig.getParameter<double>("coneSize");
0087
0088 mExtrapolations = iConfig.getParameter<edm::InputTag>("extrapolations");
0089
0090 mJPTalgo = new JetPlusTrackCorrector(iConfig, consumesCollector());
0091 if (useZSP)
0092 mZSPalgo = new ZSPJPTJetCorrector(iConfig);
0093
0094 produces<reco::JPTJetCollection>().setBranchAlias(alias);
0095
0096 input_jets_token_ = consumes<edm::View<reco::CaloJet> >(src);
0097 input_vertex_token_ = consumes<reco::VertexCollection>(srcPVs_);
0098 input_tracks_token_ = consumes<reco::TrackCollection>(mTracks);
0099 input_extrapolations_token_ = consumes<std::vector<reco::TrackExtrapolation> >(mExtrapolations);
0100 }
0101
0102 JetPlusTrackProducerAA::~JetPlusTrackProducerAA() {
0103
0104
0105 }
0106
0107
0108
0109
0110
0111
0112 void JetPlusTrackProducerAA::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0113 using namespace edm;
0114
0115
0116
0117 edm::Handle<reco::TrackCollection> tracks_h;
0118 iEvent.getByToken(input_tracks_token_, tracks_h);
0119
0120 auto const& jets_h = iEvent.get(input_jets_token_);
0121
0122 std::vector<reco::TrackRef> fTracks;
0123 fTracks.reserve(tracks_h->size());
0124 for (unsigned i = 0; i < tracks_h->size(); ++i) {
0125 fTracks.push_back(reco::TrackRef(tracks_h, i));
0126 }
0127
0128 edm::Handle<std::vector<reco::TrackExtrapolation> > extrapolations_h;
0129 iEvent.getByToken(input_extrapolations_token_, extrapolations_h);
0130
0131 auto pOut = std::make_unique<reco::JPTJetCollection>();
0132
0133 reco::JPTJetCollection tmpColl;
0134
0135 int iJet = 0;
0136 for (auto const& oldjet : jets_h) {
0137 reco::CaloJet corrected = oldjet;
0138
0139
0140
0141 double factorZSP = 1.;
0142 if (useZSP)
0143 factorZSP = mZSPalgo->correction(corrected, iEvent, iSetup);
0144
0145 corrected.scaleEnergy(factorZSP);
0146
0147
0148
0149 double scaleJPT = 1.;
0150
0151 math::XYZTLorentzVector p4;
0152
0153
0154 jpt::MatchedTracks pions;
0155 jpt::MatchedTracks muons;
0156 jpt::MatchedTracks elecs;
0157 bool validMatches = false;
0158
0159 if (!vectorial_) {
0160 scaleJPT = mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, pions, muons, elecs, validMatches);
0161 p4 = math::XYZTLorentzVector(corrected.px() * scaleJPT,
0162 corrected.py() * scaleJPT,
0163 corrected.pz() * scaleJPT,
0164 corrected.energy() * scaleJPT);
0165 } else {
0166 scaleJPT = mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, p4, pions, muons, elecs, validMatches);
0167 }
0168
0169 reco::JPTJet::Specific specific;
0170
0171 if (validMatches) {
0172 specific.pionsInVertexInCalo = pions.inVertexInCalo_;
0173 specific.pionsInVertexOutCalo = pions.inVertexOutOfCalo_;
0174 specific.pionsOutVertexInCalo = pions.outOfVertexInCalo_;
0175 specific.muonsInVertexInCalo = muons.inVertexInCalo_;
0176 specific.muonsInVertexOutCalo = muons.inVertexOutOfCalo_;
0177 specific.muonsOutVertexInCalo = muons.outOfVertexInCalo_;
0178 specific.elecsInVertexInCalo = elecs.inVertexInCalo_;
0179 specific.elecsInVertexOutCalo = elecs.inVertexOutOfCalo_;
0180 specific.elecsOutVertexInCalo = elecs.outOfVertexInCalo_;
0181 }
0182
0183
0184 specific.theCaloJetRef = edm::RefToBase<reco::Jet>(jets_h.refAt(iJet));
0185 specific.mResponseOfChargedWithEff = (float)mJPTalgo->getResponseOfChargedWithEff();
0186 specific.mResponseOfChargedWithoutEff = (float)mJPTalgo->getResponseOfChargedWithoutEff();
0187 specific.mSumPtOfChargedWithEff = (float)mJPTalgo->getSumPtWithEff();
0188 specific.mSumPtOfChargedWithoutEff = (float)mJPTalgo->getSumPtWithoutEff();
0189 specific.mSumEnergyOfChargedWithEff = (float)mJPTalgo->getSumEnergyWithEff();
0190 specific.mSumEnergyOfChargedWithoutEff = (float)mJPTalgo->getSumEnergyWithoutEff();
0191 specific.mChargedHadronEnergy = (float)mJPTalgo->getSumEnergyWithoutEff();
0192
0193
0194 double deR2Tr = 0.;
0195 double deEta2Tr = 0.;
0196 double dePhi2Tr = 0.;
0197 double Zch = 0.;
0198 double Pout2 = 0.;
0199 double Pout = 0.;
0200 double denominator_tracks = 0.;
0201 int ntracks = 0;
0202
0203 for (reco::TrackRefVector::const_iterator it = pions.inVertexInCalo_.begin(); it != pions.inVertexInCalo_.end();
0204 it++) {
0205 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0206 double deEta = (*it)->eta() - p4.eta();
0207 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0208 if ((**it).ptError() / (**it).pt() < 0.1) {
0209 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0210 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0211 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0212 denominator_tracks = denominator_tracks + (*it)->pt();
0213 Zch = Zch + (*it)->pt();
0214
0215 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0216 ntracks++;
0217 }
0218 }
0219 for (reco::TrackRefVector::const_iterator it = muons.inVertexInCalo_.begin(); it != muons.inVertexInCalo_.end();
0220 it++) {
0221 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0222 double deEta = (*it)->eta() - p4.eta();
0223 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0224 if ((**it).ptError() / (**it).pt() < 0.1) {
0225 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0226 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0227 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0228 denominator_tracks = denominator_tracks + (*it)->pt();
0229 Zch = Zch + (*it)->pt();
0230
0231 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0232 ntracks++;
0233 }
0234 }
0235 for (reco::TrackRefVector::const_iterator it = elecs.inVertexInCalo_.begin(); it != elecs.inVertexInCalo_.end();
0236 it++) {
0237 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0238 double deEta = (*it)->eta() - p4.eta();
0239 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0240 if ((**it).ptError() / (**it).pt() < 0.1) {
0241 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0242 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0243 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0244 denominator_tracks = denominator_tracks + (*it)->pt();
0245 Zch = Zch + (*it)->pt();
0246
0247 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0248 ntracks++;
0249 }
0250 }
0251 for (reco::TrackRefVector::const_iterator it = pions.inVertexOutOfCalo_.begin();
0252 it != pions.inVertexOutOfCalo_.end();
0253 it++) {
0254 Zch = Zch + (*it)->pt();
0255 }
0256 for (reco::TrackRefVector::const_iterator it = muons.inVertexOutOfCalo_.begin();
0257 it != muons.inVertexOutOfCalo_.end();
0258 it++) {
0259 Zch = Zch + (*it)->pt();
0260 }
0261 for (reco::TrackRefVector::const_iterator it = elecs.inVertexOutOfCalo_.begin();
0262 it != elecs.inVertexOutOfCalo_.end();
0263 it++) {
0264 Zch = Zch + (*it)->pt();
0265 }
0266
0267 if (mJPTalgo->getSumPtForBeta() > 0.)
0268 Zch = Zch / mJPTalgo->getSumPtForBeta();
0269
0270 if (ntracks > 0) {
0271 Pout = sqrt(fabs(Pout2)) / ntracks;
0272 }
0273
0274 if (denominator_tracks != 0) {
0275 deR2Tr = deR2Tr / denominator_tracks;
0276 deEta2Tr = deEta2Tr / denominator_tracks;
0277 dePhi2Tr = dePhi2Tr / denominator_tracks;
0278 }
0279
0280 specific.R2momtr = deR2Tr;
0281 specific.Eta2momtr = deEta2Tr;
0282 specific.Phi2momtr = dePhi2Tr;
0283 specific.Pout = Pout;
0284 specific.Zch = Zch;
0285
0286
0287 reco::Particle::Point vertex_ = reco::Jet::Point(0, 0, 0);
0288
0289
0290 edm::Handle<reco::VertexCollection> pvCollection;
0291 iEvent.getByToken(input_vertex_token_, pvCollection);
0292 if (pvCollection.isValid() && !pvCollection->empty())
0293 vertex_ = pvCollection->begin()->position();
0294
0295 reco::JPTJet fJet(p4, vertex_, specific, corrected.getJetConstituents());
0296
0297
0298
0299 iJet++;
0300 tmpColl.push_back(fJet);
0301 }
0302
0303
0304
0305
0306 reco::TrackRefVector trBgOutOfCalo;
0307 reco::TrackRefVector trBgOutOfVertex = calculateBGtracksJet(tmpColl, fTracks, extrapolations_h, trBgOutOfCalo);
0308
0309
0310 std::map<reco::JPTJetCollection::iterator, double> AreaNonJet;
0311
0312 for (reco::JPTJetCollection::iterator ij1 = tmpColl.begin(); ij1 != tmpColl.end(); ij1++) {
0313 int nj1 = 1;
0314 for (reco::JPTJetCollection::iterator ij2 = tmpColl.begin(); ij2 != tmpColl.end(); ij2++) {
0315 if (ij2 == ij1)
0316 continue;
0317 if (fabs((*ij1).eta() - (*ij2).eta()) > 0.5)
0318 continue;
0319 nj1++;
0320 }
0321
0322 AreaNonJet[ij1] = 4 * M_PI * mConeSize - nj1 * 4 * mConeSize * mConeSize;
0323 }
0324
0325
0326
0327 for (reco::JPTJetCollection::iterator ij = tmpColl.begin(); ij != tmpColl.end(); ij++) {
0328
0329
0330 const reco::TrackRefVector pioninin = (*ij).getPionsInVertexInCalo();
0331 const reco::TrackRefVector pioninout = (*ij).getPionsInVertexOutCalo();
0332
0333 double ja = (AreaNonJet.find(ij))->second;
0334
0335 double factorPU = mJPTalgo->correctAA(*ij, trBgOutOfVertex, mConeSize, pioninin, pioninout, ja, trBgOutOfCalo);
0336
0337 (*ij).scaleEnergy(factorPU);
0338
0339
0340 pOut->push_back(*ij);
0341 }
0342
0343 iEvent.put(std::move(pOut));
0344 }
0345
0346
0347
0348
0349 reco::TrackRefVector JetPlusTrackProducerAA::calculateBGtracksJet(
0350 reco::JPTJetCollection& fJets,
0351 std::vector<reco::TrackRef>& fTracks,
0352 edm::Handle<std::vector<reco::TrackExtrapolation> >& extrapolations_h,
0353 reco::TrackRefVector& trBgOutOfCalo) {
0354 reco::TrackRefVector trBgOutOfVertex;
0355
0356 for (unsigned t = 0; t < fTracks.size(); ++t) {
0357 int track_bg = 0;
0358
0359 const reco::Track* track = &*(fTracks[t]);
0360 double trackEta = track->eta();
0361 double trackPhi = track->phi();
0362
0363
0364 for (unsigned j = 0; j < fJets.size(); ++j) {
0365 const reco::Jet* jet = &(fJets[j]);
0366 double jetEta = jet->eta();
0367 double jetPhi = jet->phi();
0368
0369 if (fabs(jetEta - trackEta) < mConeSize) {
0370 double dphiTrackJet = deltaPhi(trackPhi, jetPhi);
0371 if (dphiTrackJet < mConeSize) {
0372 track_bg = 1;
0373 }
0374 }
0375 }
0376
0377 if (track_bg == 0) {
0378 trBgOutOfVertex.push_back(fTracks[t]);
0379 }
0380
0381 }
0382
0383
0384 int nValid = 0;
0385 for (std::vector<reco::TrackExtrapolation>::const_iterator xtrpBegin = extrapolations_h->begin(),
0386 xtrpEnd = extrapolations_h->end(),
0387 ixtrp = xtrpBegin;
0388 ixtrp != xtrpEnd;
0389 ++ixtrp) {
0390 nValid++;
0391
0392 reco::TrackRefVector::iterator it = find(trBgOutOfVertex.begin(), trBgOutOfVertex.end(), (*ixtrp).track());
0393
0394 if (it != trBgOutOfVertex.end()) {
0395 trBgOutOfCalo.push_back(*it);
0396 }
0397 }
0398
0399 return trBgOutOfVertex;
0400 }
0401
0402
0403
