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 #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 math::XYZTLorentzVector p4;
0148
0149
0150 jpt::MatchedTracks pions;
0151 jpt::MatchedTracks muons;
0152 jpt::MatchedTracks elecs;
0153 bool validMatches = false;
0154
0155 if (!vectorial_) {
0156 double scaleJPT = mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, pions, muons, elecs, validMatches);
0157 p4 = math::XYZTLorentzVector(corrected.px() * scaleJPT,
0158 corrected.py() * scaleJPT,
0159 corrected.pz() * scaleJPT,
0160 corrected.energy() * scaleJPT);
0161 } else {
0162 mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, p4, pions, muons, elecs, validMatches);
0163 }
0164
0165 reco::JPTJet::Specific specific;
0166
0167 if (validMatches) {
0168 specific.pionsInVertexInCalo = pions.inVertexInCalo_;
0169 specific.pionsInVertexOutCalo = pions.inVertexOutOfCalo_;
0170 specific.pionsOutVertexInCalo = pions.outOfVertexInCalo_;
0171 specific.muonsInVertexInCalo = muons.inVertexInCalo_;
0172 specific.muonsInVertexOutCalo = muons.inVertexOutOfCalo_;
0173 specific.muonsOutVertexInCalo = muons.outOfVertexInCalo_;
0174 specific.elecsInVertexInCalo = elecs.inVertexInCalo_;
0175 specific.elecsInVertexOutCalo = elecs.inVertexOutOfCalo_;
0176 specific.elecsOutVertexInCalo = elecs.outOfVertexInCalo_;
0177 }
0178
0179
0180 specific.theCaloJetRef = edm::RefToBase<reco::Jet>(jets_h.refAt(iJet));
0181 specific.mResponseOfChargedWithEff = (float)mJPTalgo->getResponseOfChargedWithEff();
0182 specific.mResponseOfChargedWithoutEff = (float)mJPTalgo->getResponseOfChargedWithoutEff();
0183 specific.mSumPtOfChargedWithEff = (float)mJPTalgo->getSumPtWithEff();
0184 specific.mSumPtOfChargedWithoutEff = (float)mJPTalgo->getSumPtWithoutEff();
0185 specific.mSumEnergyOfChargedWithEff = (float)mJPTalgo->getSumEnergyWithEff();
0186 specific.mSumEnergyOfChargedWithoutEff = (float)mJPTalgo->getSumEnergyWithoutEff();
0187 specific.mChargedHadronEnergy = (float)mJPTalgo->getSumEnergyWithoutEff();
0188
0189
0190 double deR2Tr = 0.;
0191 double deEta2Tr = 0.;
0192 double dePhi2Tr = 0.;
0193 double Zch = 0.;
0194 double Pout2 = 0.;
0195 double Pout = 0.;
0196 double denominator_tracks = 0.;
0197 int ntracks = 0;
0198
0199 for (reco::TrackRefVector::const_iterator it = pions.inVertexInCalo_.begin(); it != pions.inVertexInCalo_.end();
0200 it++) {
0201 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0202 double deEta = (*it)->eta() - p4.eta();
0203 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0204 if ((**it).ptError() / (**it).pt() < 0.1) {
0205 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0206 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0207 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0208 denominator_tracks = denominator_tracks + (*it)->pt();
0209 Zch = Zch + (*it)->pt();
0210
0211 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0212 ntracks++;
0213 }
0214 }
0215 for (reco::TrackRefVector::const_iterator it = muons.inVertexInCalo_.begin(); it != muons.inVertexInCalo_.end();
0216 it++) {
0217 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0218 double deEta = (*it)->eta() - p4.eta();
0219 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0220 if ((**it).ptError() / (**it).pt() < 0.1) {
0221 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0222 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0223 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0224 denominator_tracks = denominator_tracks + (*it)->pt();
0225 Zch = Zch + (*it)->pt();
0226
0227 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0228 ntracks++;
0229 }
0230 }
0231 for (reco::TrackRefVector::const_iterator it = elecs.inVertexInCalo_.begin(); it != elecs.inVertexInCalo_.end();
0232 it++) {
0233 double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0234 double deEta = (*it)->eta() - p4.eta();
0235 double dePhi = deltaPhi((*it)->phi(), p4.phi());
0236 if ((**it).ptError() / (**it).pt() < 0.1) {
0237 deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0238 deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0239 dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0240 denominator_tracks = denominator_tracks + (*it)->pt();
0241 Zch = Zch + (*it)->pt();
0242
0243 Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0244 ntracks++;
0245 }
0246 }
0247 for (reco::TrackRefVector::const_iterator it = pions.inVertexOutOfCalo_.begin();
0248 it != pions.inVertexOutOfCalo_.end();
0249 it++) {
0250 Zch = Zch + (*it)->pt();
0251 }
0252 for (reco::TrackRefVector::const_iterator it = muons.inVertexOutOfCalo_.begin();
0253 it != muons.inVertexOutOfCalo_.end();
0254 it++) {
0255 Zch = Zch + (*it)->pt();
0256 }
0257 for (reco::TrackRefVector::const_iterator it = elecs.inVertexOutOfCalo_.begin();
0258 it != elecs.inVertexOutOfCalo_.end();
0259 it++) {
0260 Zch = Zch + (*it)->pt();
0261 }
0262
0263 if (mJPTalgo->getSumPtForBeta() > 0.)
0264 Zch = Zch / mJPTalgo->getSumPtForBeta();
0265
0266 if (ntracks > 0) {
0267 Pout = sqrt(fabs(Pout2)) / ntracks;
0268 }
0269
0270 if (denominator_tracks != 0) {
0271 deR2Tr = deR2Tr / denominator_tracks;
0272 deEta2Tr = deEta2Tr / denominator_tracks;
0273 dePhi2Tr = dePhi2Tr / denominator_tracks;
0274 }
0275
0276 specific.R2momtr = deR2Tr;
0277 specific.Eta2momtr = deEta2Tr;
0278 specific.Phi2momtr = dePhi2Tr;
0279 specific.Pout = Pout;
0280 specific.Zch = Zch;
0281
0282
0283 reco::Particle::Point vertex_ = reco::Jet::Point(0, 0, 0);
0284
0285
0286 edm::Handle<reco::VertexCollection> pvCollection;
0287 iEvent.getByToken(input_vertex_token_, pvCollection);
0288 if (pvCollection.isValid() && !pvCollection->empty())
0289 vertex_ = pvCollection->begin()->position();
0290
0291 reco::JPTJet fJet(p4, vertex_, specific, corrected.getJetConstituents());
0292
0293
0294
0295 iJet++;
0296 tmpColl.push_back(fJet);
0297 }
0298
0299
0300
0301
0302 reco::TrackRefVector trBgOutOfCalo;
0303 reco::TrackRefVector trBgOutOfVertex = calculateBGtracksJet(tmpColl, fTracks, extrapolations_h, trBgOutOfCalo);
0304
0305
0306 std::map<reco::JPTJetCollection::iterator, double> AreaNonJet;
0307
0308 for (reco::JPTJetCollection::iterator ij1 = tmpColl.begin(); ij1 != tmpColl.end(); ij1++) {
0309 int nj1 = 1;
0310 for (reco::JPTJetCollection::iterator ij2 = tmpColl.begin(); ij2 != tmpColl.end(); ij2++) {
0311 if (ij2 == ij1)
0312 continue;
0313 if (fabs((*ij1).eta() - (*ij2).eta()) > 0.5)
0314 continue;
0315 nj1++;
0316 }
0317
0318 AreaNonJet[ij1] = 4 * M_PI * mConeSize - nj1 * 4 * mConeSize * mConeSize;
0319 }
0320
0321
0322
0323 for (reco::JPTJetCollection::iterator ij = tmpColl.begin(); ij != tmpColl.end(); ij++) {
0324
0325
0326 const reco::TrackRefVector pioninin = (*ij).getPionsInVertexInCalo();
0327 const reco::TrackRefVector pioninout = (*ij).getPionsInVertexOutCalo();
0328
0329 double ja = (AreaNonJet.find(ij))->second;
0330
0331 double factorPU = mJPTalgo->correctAA(*ij, trBgOutOfVertex, mConeSize, pioninin, pioninout, ja, trBgOutOfCalo);
0332
0333 (*ij).scaleEnergy(factorPU);
0334
0335
0336 pOut->push_back(*ij);
0337 }
0338
0339 iEvent.put(std::move(pOut));
0340 }
0341
0342
0343
0344
0345 reco::TrackRefVector JetPlusTrackProducerAA::calculateBGtracksJet(
0346 reco::JPTJetCollection& fJets,
0347 std::vector<reco::TrackRef>& fTracks,
0348 edm::Handle<std::vector<reco::TrackExtrapolation> >& extrapolations_h,
0349 reco::TrackRefVector& trBgOutOfCalo) {
0350 reco::TrackRefVector trBgOutOfVertex;
0351
0352 for (unsigned t = 0; t < fTracks.size(); ++t) {
0353 int track_bg = 0;
0354
0355 const reco::Track* track = &*(fTracks[t]);
0356 double trackEta = track->eta();
0357 double trackPhi = track->phi();
0358
0359
0360 for (unsigned j = 0; j < fJets.size(); ++j) {
0361 const reco::Jet* jet = &(fJets[j]);
0362 double jetEta = jet->eta();
0363 double jetPhi = jet->phi();
0364
0365 if (fabs(jetEta - trackEta) < mConeSize) {
0366 double dphiTrackJet = deltaPhi(trackPhi, jetPhi);
0367 if (dphiTrackJet < mConeSize) {
0368 track_bg = 1;
0369 }
0370 }
0371 }
0372
0373 if (track_bg == 0) {
0374 trBgOutOfVertex.push_back(fTracks[t]);
0375 }
0376
0377 }
0378
0379
0380 for (std::vector<reco::TrackExtrapolation>::const_iterator xtrpBegin = extrapolations_h->begin(),
0381 xtrpEnd = extrapolations_h->end(),
0382 ixtrp = xtrpBegin;
0383 ixtrp != xtrpEnd;
0384 ++ixtrp) {
0385 reco::TrackRefVector::iterator it = find(trBgOutOfVertex.begin(), trBgOutOfVertex.end(), (*ixtrp).track());
0386
0387 if (it != trBgOutOfVertex.end()) {
0388 trBgOutOfCalo.push_back(*it);
0389 }
0390 }
0391
0392 return trBgOutOfVertex;
0393 }
0394
0395
0396
