Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-10-08 23:09:56

0001 // -*- C++ -*-
0002 //
0003 // Package:    JetPlusTrack
0004 // Class:      JetPlusTrack
0005 //
0006 /**\class JetPlusTrackProducer JetPlusTrackProducer.cc JetPlusTrackProducer.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Olga Kodolova,40 R-A12,+41227671273,
0015 //         Created:  Fri Feb 19 10:14:02 CET 2010
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
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 // constants, enums and typedefs
0067 //
0068 
0069 //
0070 // static data member definitions
0071 //
0072 
0073 //
0074 // constructors and destructor
0075 //
0076 JetPlusTrackProducerAA::JetPlusTrackProducerAA(const edm::ParameterSet& iConfig) {
0077   //register your products
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   // do anything here that needs to be done at desctruction time
0104   // (e.g. close files, deallocate resources etc.)
0105 }
0106 
0107 //
0108 // member functions
0109 //
0110 
0111 // ------------ method called to produce the data  ------------
0112 void JetPlusTrackProducerAA::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0113   using namespace edm;
0114 
0115   // get stuff from Event
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     // ZSP corrections
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     // Construct JPTJet constituent
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     // Fill JPT Specific
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     // Fill Charged Jet shape parameters
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     // Create JPT jet
0283     reco::Particle::Point vertex_ = reco::Jet::Point(0, 0, 0);
0284 
0285     // If we add primary vertex
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     // Temporarily collection before correction for background
0294 
0295     iJet++;
0296     tmpColl.push_back(fJet);
0297   }
0298 
0299   //=======================================================================================================>
0300   // Correction for background
0301 
0302   reco::TrackRefVector trBgOutOfCalo;
0303   reco::TrackRefVector trBgOutOfVertex = calculateBGtracksJet(tmpColl, fTracks, extrapolations_h, trBgOutOfCalo);
0304 
0305   //===> Area without Jets
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     // Correct JPTjet for background tracks
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     // Output module
0336     pOut->push_back(*ij);
0337   }
0338 
0339   iEvent.put(std::move(pOut));
0340 }
0341 // -----------------------------------------------
0342 // ------------ calculateBGtracksJet  ------------
0343 // ------------ Tracks not included in jets ------
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     //loop on jets
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     }  //jets
0372 
0373     if (track_bg == 0) {
0374       trBgOutOfVertex.push_back(fTracks[t]);
0375     }
0376 
0377   }  //tracks
0378 
0379   //=====> Propagate BG tracks to calo
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 //define this as a plug-in
0396 //DEFINE_FWK_MODULE(JetPlusTrackProducerAA);