Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:59:48

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     // JPT corrections
0148 
0149     double scaleJPT = 1.;
0150 
0151     math::XYZTLorentzVector p4;
0152 
0153     // Construct JPTJet constituent
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     // Fill JPT Specific
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     // Fill Charged Jet shape parameters
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     // Create JPT jet
0287     reco::Particle::Point vertex_ = reco::Jet::Point(0, 0, 0);
0288 
0289     // If we add primary vertex
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     // Temporarily collection before correction for background
0298 
0299     iJet++;
0300     tmpColl.push_back(fJet);
0301   }
0302 
0303   //=======================================================================================================>
0304   // Correction for background
0305 
0306   reco::TrackRefVector trBgOutOfCalo;
0307   reco::TrackRefVector trBgOutOfVertex = calculateBGtracksJet(tmpColl, fTracks, extrapolations_h, trBgOutOfCalo);
0308 
0309   //===> Area without Jets
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     // Correct JPTjet for background tracks
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     // Output module
0340     pOut->push_back(*ij);
0341   }
0342 
0343   iEvent.put(std::move(pOut));
0344 }
0345 // -----------------------------------------------
0346 // ------------ calculateBGtracksJet  ------------
0347 // ------------ Tracks not included in jets ------
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     //loop on jets
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     }  //jets
0376 
0377     if (track_bg == 0) {
0378       trBgOutOfVertex.push_back(fTracks[t]);
0379     }
0380 
0381   }  //tracks
0382 
0383   //=====> Propagate BG tracks to calo
0384   for (std::vector<reco::TrackExtrapolation>::const_iterator xtrpBegin = extrapolations_h->begin(),
0385                                                              xtrpEnd = extrapolations_h->end(),
0386                                                              ixtrp = xtrpBegin;
0387        ixtrp != xtrpEnd;
0388        ++ixtrp) {
0389     reco::TrackRefVector::iterator it = find(trBgOutOfVertex.begin(), trBgOutOfVertex.end(), (*ixtrp).track());
0390 
0391     if (it != trBgOutOfVertex.end()) {
0392       trBgOutOfCalo.push_back(*it);
0393     }
0394   }
0395 
0396   return trBgOutOfVertex;
0397 }
0398 
0399 //define this as a plug-in
0400 //DEFINE_FWK_MODULE(JetPlusTrackProducerAA);