Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 // -*- C++ -*-
0002 //
0003 // Package:    JetPlusTracks
0004 // Class:      JetPlusTrackProducer
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 
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 // constants, enums and typedefs
0051 //
0052 
0053 //
0054 // static data member definitions
0055 //
0056 
0057 //
0058 // constructors and destructor
0059 //
0060 JetPlusTrackProducer::JetPlusTrackProducer(const edm::ParameterSet& iConfig) {
0061   //register your products
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   // do anything here that needs to be done at desctruction time
0089   // (e.g. close files, deallocate resources etc.)
0090 }
0091 
0092 //
0093 // member functions
0094 //
0095 bool sort_by_pt(const reco::JPTJet& a, const reco::JPTJet& b) { return (a.pt() > b.pt()); }
0096 
0097 // ------------ method called to produce the data  ------------
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   double scaleJPT = 1.;
0111   for (auto const& jet : iEvent.get(input_trackjets_token_)) {
0112     int icalo = -1;
0113     int i = 0;
0114     for (auto const& oldjet : addjets_h) {
0115       double dr2 = deltaR2(jet, oldjet);
0116       if (dr2 <= dRcone_ * dRcone_) {
0117         icalo = i;
0118       }
0119       i++;
0120     }  // Calojets
0121     if (icalo < 0)
0122       continue;
0123     auto const& mycalo = addjets_h[icalo];
0124     std::vector<edm::Ptr<reco::Track> > tracksinjet = jet.tracks();
0125     reco::TrackRefVector tracksincalo;
0126     reco::TrackRefVector tracksinvert;
0127     for (auto const& itrack : tracksinjet) {
0128       for (auto const& ixtrp : iExtrapolations) {
0129         if (ixtrp.positions().empty())
0130           continue;
0131         if (usePAT_) {
0132           double mydphi = deltaPhi(ixtrp.track()->phi(), itrack->phi());
0133           if (fabs(ixtrp.track()->pt() - itrack->pt()) > 0.001 || fabs(ixtrp.track()->eta() - itrack->eta()) > 0.001 ||
0134               mydphi > 0.001)
0135             continue;
0136         } else {
0137           if (itrack.id() != ixtrp.track().id() || itrack.key() != ixtrp.track().key())
0138             continue;
0139         }
0140         tracksinvert.push_back(ixtrp.track());
0141         reco::TrackBase::Point const& point = ixtrp.positions().at(0);
0142         double dr2 = deltaR2(jet, point);
0143         if (dr2 <= dRcone_ * dRcone_) {
0144           tracksincalo.push_back(ixtrp.track());
0145         }
0146       }  // Track extrapolations
0147     }    // tracks
0148 
0149     const reco::TrackJet& corrected = jet;
0150     math::XYZTLorentzVector p4;
0151     jpt::MatchedTracks pions;
0152     jpt::MatchedTracks muons;
0153     jpt::MatchedTracks elecs;
0154 
0155     scaleJPT =
0156         mJPTalgo->correction(corrected, mycalo, iEvent, iSetup, tracksinvert, tracksincalo, p4, pions, muons, elecs);
0157     if (p4.pt() > ptCUT_) {
0158       reco::JPTJet::Specific jptspe;
0159       jptspe.pionsInVertexInCalo = pions.inVertexInCalo_;
0160       jptspe.pionsInVertexOutCalo = pions.inVertexOutOfCalo_;
0161       jptspe.pionsOutVertexInCalo = pions.outOfVertexInCalo_;
0162       jptspe.muonsInVertexInCalo = muons.inVertexInCalo_;
0163       jptspe.muonsInVertexOutCalo = muons.inVertexOutOfCalo_;
0164       jptspe.muonsOutVertexInCalo = muons.outOfVertexInCalo_;
0165       jptspe.elecsInVertexInCalo = elecs.inVertexInCalo_;
0166       jptspe.elecsInVertexOutCalo = elecs.inVertexOutOfCalo_;
0167       jptspe.elecsOutVertexInCalo = elecs.outOfVertexInCalo_;
0168       reco::CaloJetRef myjet(pOut1RefProd, idxCaloJet++);
0169       jptspe.theCaloJetRef = edm::RefToBase<reco::Jet>(myjet);
0170       jptspe.JPTSeed = 1;
0171       reco::JPTJet fJet(p4, jet.primaryVertex()->position(), jptspe, mycalo.getJetConstituents());
0172       pOut->push_back(fJet);
0173       pOut1->push_back(mycalo);
0174     }
0175   }  // trackjets
0176 
0177   int iJet = 0;
0178   for (auto const& oldjet : jets_h) {
0179     reco::CaloJet corrected = oldjet;
0180 
0181     // ZSP corrections
0182     double factorZSP = 1.;
0183     if (useZSP_)
0184       factorZSP = mZSPalgo->correction(corrected, iEvent, iSetup);
0185     corrected.scaleEnergy(factorZSP);
0186 
0187     // JPT corrections
0188     scaleJPT = 1.;
0189 
0190     math::XYZTLorentzVector p4;
0191 
0192     jpt::MatchedTracks pions;
0193     jpt::MatchedTracks muons;
0194     jpt::MatchedTracks elecs;
0195     bool validMatches = false;
0196 
0197     if (!vectorial_) {
0198       scaleJPT = mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, pions, muons, elecs, validMatches);
0199       p4 = math::XYZTLorentzVector(corrected.px() * scaleJPT,
0200                                    corrected.py() * scaleJPT,
0201                                    corrected.pz() * scaleJPT,
0202                                    corrected.energy() * scaleJPT);
0203     } else {
0204       scaleJPT = mJPTalgo->correction(corrected, oldjet, iEvent, iSetup, p4, pions, muons, elecs, validMatches);
0205     }
0206 
0207     reco::JPTJet::Specific specific;
0208 
0209     if (validMatches) {
0210       specific.pionsInVertexInCalo = pions.inVertexInCalo_;
0211       specific.pionsInVertexOutCalo = pions.inVertexOutOfCalo_;
0212       specific.pionsOutVertexInCalo = pions.outOfVertexInCalo_;
0213       specific.muonsInVertexInCalo = muons.inVertexInCalo_;
0214       specific.muonsInVertexOutCalo = muons.inVertexOutOfCalo_;
0215       specific.muonsOutVertexInCalo = muons.outOfVertexInCalo_;
0216       specific.elecsInVertexInCalo = elecs.inVertexInCalo_;
0217       specific.elecsInVertexOutCalo = elecs.inVertexOutOfCalo_;
0218       specific.elecsOutVertexInCalo = elecs.outOfVertexInCalo_;
0219     }
0220 
0221     // Fill JPT Specific
0222     specific.theCaloJetRef = edm::RefToBase<reco::Jet>(jets_h.refAt(iJet));
0223     specific.mResponseOfChargedWithEff = (float)mJPTalgo->getResponseOfChargedWithEff();
0224     specific.mResponseOfChargedWithoutEff = (float)mJPTalgo->getResponseOfChargedWithoutEff();
0225     specific.mSumPtOfChargedWithEff = (float)mJPTalgo->getSumPtWithEff();
0226     specific.mSumPtOfChargedWithoutEff = (float)mJPTalgo->getSumPtWithoutEff();
0227     specific.mSumEnergyOfChargedWithEff = (float)mJPTalgo->getSumEnergyWithEff();
0228     specific.mSumEnergyOfChargedWithoutEff = (float)mJPTalgo->getSumEnergyWithoutEff();
0229     specific.mChargedHadronEnergy = (float)mJPTalgo->getSumEnergyWithoutEff();
0230 
0231     // Fill Charged Jet shape parameters
0232 
0233     double deR2Tr = 0.;
0234     double deEta2Tr = 0.;
0235     double dePhi2Tr = 0.;
0236     double Zch = 0.;
0237     double Pout2 = 0.;
0238     double Pout = 0.;
0239     double denominator_tracks = 0.;
0240     int ntracks = 0;
0241 
0242     for (reco::TrackRefVector::const_iterator it = pions.inVertexInCalo_.begin(); it != pions.inVertexInCalo_.end();
0243          it++) {
0244       double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0245       double deEta = (*it)->eta() - p4.eta();
0246       double dePhi = deltaPhi((*it)->phi(), p4.phi());
0247       if ((**it).ptError() / (**it).pt() < 0.1) {
0248         deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0249         deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0250         dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0251         denominator_tracks = denominator_tracks + (*it)->pt();
0252         Zch = Zch + (*it)->pt();
0253 
0254         Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0255         ntracks++;
0256       }
0257     }
0258     for (reco::TrackRefVector::const_iterator it = muons.inVertexInCalo_.begin(); it != muons.inVertexInCalo_.end();
0259          it++) {
0260       double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0261       double deEta = (*it)->eta() - p4.eta();
0262       double dePhi = deltaPhi((*it)->phi(), p4.phi());
0263       if ((**it).ptError() / (**it).pt() < 0.1) {
0264         deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0265         deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0266         dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0267         denominator_tracks = denominator_tracks + (*it)->pt();
0268         Zch = Zch + (*it)->pt();
0269 
0270         Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0271         ntracks++;
0272       }
0273     }
0274     for (reco::TrackRefVector::const_iterator it = elecs.inVertexInCalo_.begin(); it != elecs.inVertexInCalo_.end();
0275          it++) {
0276       double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
0277       double deEta = (*it)->eta() - p4.eta();
0278       double dePhi = deltaPhi((*it)->phi(), p4.phi());
0279       if ((**it).ptError() / (**it).pt() < 0.1) {
0280         deR2Tr = deR2Tr + deR * deR * (*it)->pt();
0281         deEta2Tr = deEta2Tr + deEta * deEta * (*it)->pt();
0282         dePhi2Tr = dePhi2Tr + dePhi * dePhi * (*it)->pt();
0283         denominator_tracks = denominator_tracks + (*it)->pt();
0284         Zch = Zch + (*it)->pt();
0285 
0286         Pout2 = Pout2 + (**it).p() * (**it).p() - (Zch * p4.P()) * (Zch * p4.P());
0287         ntracks++;
0288       }
0289     }
0290 
0291     for (reco::TrackRefVector::const_iterator it = pions.inVertexOutOfCalo_.begin();
0292          it != pions.inVertexOutOfCalo_.end();
0293          it++) {
0294       Zch = Zch + (*it)->pt();
0295     }
0296     for (reco::TrackRefVector::const_iterator it = muons.inVertexOutOfCalo_.begin();
0297          it != muons.inVertexOutOfCalo_.end();
0298          it++) {
0299       Zch = Zch + (*it)->pt();
0300     }
0301     for (reco::TrackRefVector::const_iterator it = elecs.inVertexOutOfCalo_.begin();
0302          it != elecs.inVertexOutOfCalo_.end();
0303          it++) {
0304       Zch = Zch + (*it)->pt();
0305     }
0306 
0307     if (mJPTalgo->getSumPtForBeta() > 0.)
0308       Zch = Zch / mJPTalgo->getSumPtForBeta();
0309 
0310     if (ntracks > 0) {
0311       Pout = sqrt(fabs(Pout2)) / ntracks;
0312     }
0313     if (denominator_tracks != 0) {
0314       deR2Tr = deR2Tr / denominator_tracks;
0315       deEta2Tr = deEta2Tr / denominator_tracks;
0316       dePhi2Tr = dePhi2Tr / denominator_tracks;
0317     }
0318 
0319     specific.R2momtr = deR2Tr;
0320     specific.Eta2momtr = deEta2Tr;
0321     specific.Phi2momtr = dePhi2Tr;
0322     specific.Pout = Pout;
0323     specific.Zch = Zch;
0324 
0325     // Create JPT jet
0326     reco::Particle::Point vertex_ = reco::Jet::Point(0, 0, 0);
0327 
0328     // If we add primary vertex
0329     edm::Handle<reco::VertexCollection> pvCollection;
0330     iEvent.getByToken(input_vertex_token_, pvCollection);
0331     if (pvCollection.isValid() && !pvCollection->empty())
0332       vertex_ = pvCollection->begin()->position();
0333 
0334     reco::JPTJet fJet(p4, vertex_, specific, corrected.getJetConstituents());
0335     iJet++;
0336 
0337     // Output module
0338     if (fJet.pt() > ptCUT_)
0339       pOut->push_back(fJet);
0340   }
0341   std::sort(pOut->begin(), pOut->end(), sort_by_pt);
0342   iEvent.put(std::move(pOut1));
0343   iEvent.put(std::move(pOut));
0344 }
0345 
0346 //define this as a plug-in
0347 //DEFINE_FWK_MODULE(JetPlusTrackProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!