Back to home page

Project CMSSW displayed by LXR
 
 

    

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

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