Skip to content

Commit

Permalink
Merge pull request #14269 from silviodonato/pixelPuJetId
Browse files Browse the repository at this point in the history 
Pixel PU jet update
  • Loading branch information
@cmsbuild
cmsbuild committed Apr 28, 2016
2 parents 5e996cd + 81d41ba commit 5fb74b0
Showing 1 changed file with 161 additions and 143 deletions.
304 changes: 161 additions & 143 deletions HLTrigger/JetMET/plugins/PixelJetPuId.cc
View file
Expand Up @@ -5,13 +5,13 @@
//
/**\class PixelJetPuId PixelJetPuId.cc RecoBTag/PixelJetPuId/src/PixelJetPuId.cc
Description:
The PixelJetPuId module select all the pixel tracks compatible with a jet.
If the sum of the tracks momentum is under a threshold the jet is tagged as "PUjets".
Description:
The PixelJetPuId module select all the pixel tracks compatible with a jet.
If the sum of the tracks momentum is under a threshold the jet is tagged as "PUjets".
Implementation:
[Notes on implementation]
*/
Implementation:
[Notes on implementation]
*/
//
// Original Author: Silvio DONATO
// Created: Wed Dec 18 10:05:40 CET 2013
Expand All @@ -28,6 +28,7 @@

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/makeRefToBaseProdFrom.h"

#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
Expand All @@ -46,40 +47,41 @@
#include "TrackingTools/Records/interface/TransientTrackRecord.h"

#include "DataFormats/Math/interface/deltaR.h"

#include "DataFormats/BTauReco/interface/JetTag.h"
#include "DataFormats/Common/interface/RefToBase.h"
//
// class declaration
//

class PixelJetPuId : public edm::global::EDProducer <>{
public:
PixelJetPuId(const edm::ParameterSet&);
virtual ~PixelJetPuId();

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
virtual void produce(edm::StreamID sid, edm::Event&, const edm::EventSetup&) const override;


// ----------member data ---------------------------
edm::InputTag m_primaryVertex;
edm::InputTag m_tracks;
edm::InputTag m_jets;
edm::EDGetTokenT<std::vector<reco::Track> > tracksToken;
edm::EDGetTokenT<edm::View<reco::CaloJet> > jetsToken;
edm::EDGetTokenT<reco::VertexCollection> primaryVertexToken;

double m_MinTrackPt;
double m_MaxTrackChi2;
double m_MaxTrackDistanceToJet;

bool m_fwjets;
double m_mineta_fwjets;
double m_minet_fwjets;

double m_MinGoodJetTrackPt;
double m_MinGoodJetTrackPtRatio;
public:
PixelJetPuId(const edm::ParameterSet&);
virtual ~PixelJetPuId();
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
private:
virtual void produce(edm::StreamID sid, edm::Event&, const edm::EventSetup&) const override;
// ----------member data ---------------------------
edm::InputTag m_primaryVertex;
edm::InputTag m_tracks;
edm::InputTag m_jets;
edm::EDGetTokenT<std::vector<reco::Track> > tracksToken;
edm::EDGetTokenT<edm::View<reco::CaloJet> > jetsToken;
edm::EDGetTokenT<edm::View<reco::Jet> > generaljetsToken;
edm::EDGetTokenT<reco::VertexCollection> primaryVertexToken;

double m_MinTrackPt;
double m_MaxTrackChi2;
double m_MaxTrackDistanceToJet;

bool m_fwjets;
double m_mineta_fwjets;
double m_minet_fwjets;

double m_MinGoodJetTrackPt;
double m_MinGoodJetTrackPtRatio;
};


Expand All @@ -88,29 +90,31 @@ class PixelJetPuId : public edm::global::EDProducer <>{
//
PixelJetPuId::PixelJetPuId(const edm::ParameterSet& iConfig)
{
//InputTag
m_tracks = iConfig.getParameter<edm::InputTag>("tracks");
tracksToken = consumes<std::vector<reco::Track> >(m_tracks);
m_jets = iConfig.getParameter<edm::InputTag>("jets");
jetsToken = consumes<edm::View<reco::CaloJet> >(m_jets);
m_primaryVertex = iConfig.getParameter<edm::InputTag>("primaryVertex");
primaryVertexToken = consumes<reco::VertexCollection>(m_primaryVertex);

//Tracks Selection
m_MinTrackPt = iConfig.getParameter<double>("MinTrackPt");
m_MaxTrackDistanceToJet = iConfig.getParameter<double>("MaxTrackDistanceToJet");
m_MaxTrackChi2 = iConfig.getParameter<double>("MaxTrackChi2");

//A jet is defined as a signal jet if Sum(trackPt) > minPt or Sum(comp.trackPt)/CaloJetPt > minPtRatio
m_MinGoodJetTrackPt = iConfig.getParameter<double>("MinGoodJetTrackPt");
m_MinGoodJetTrackPtRatio = iConfig.getParameter<double>("MinGoodJetTrackPtRatio");

m_fwjets = iConfig.getParameter<bool>("UseForwardJetsAsNoPU");
m_mineta_fwjets = iConfig.getParameter<double>("MinEtaForwardJets");
m_minet_fwjets = iConfig.getParameter<double>("MinEtForwardJets");

produces<std::vector<reco::CaloJet> >();
produces<std::vector<reco::CaloJet> >("PUjets");
//InputTag
m_tracks = iConfig.getParameter<edm::InputTag>("tracks");
tracksToken = consumes<std::vector<reco::Track> >(m_tracks);
m_jets = iConfig.getParameter<edm::InputTag>("jets");
jetsToken = consumes<edm::View<reco::CaloJet> >(m_jets);
generaljetsToken = consumes<edm::View<reco::Jet> >(m_jets);
m_primaryVertex = iConfig.getParameter<edm::InputTag>("primaryVertex");
primaryVertexToken = consumes<reco::VertexCollection>(m_primaryVertex);

//Tracks Selection
m_MinTrackPt = iConfig.getParameter<double>("MinTrackPt");
m_MaxTrackDistanceToJet = iConfig.getParameter<double>("MaxTrackDistanceToJet");
m_MaxTrackChi2 = iConfig.getParameter<double>("MaxTrackChi2");

//A jet is defined as a signal jet if Sum(trackPt) > minPt or Sum(comp.trackPt)/CaloJetPt > minPtRatio
m_MinGoodJetTrackPt = iConfig.getParameter<double>("MinGoodJetTrackPt");
m_MinGoodJetTrackPtRatio = iConfig.getParameter<double>("MinGoodJetTrackPtRatio");

m_fwjets = iConfig.getParameter<bool>("UseForwardJetsAsNoPU");
m_mineta_fwjets = iConfig.getParameter<double>("MinEtaForwardJets");
m_minet_fwjets = iConfig.getParameter<double>("MinEtForwardJets");

produces<std::vector<reco::CaloJet> >();
produces<std::vector<reco::CaloJet> >("PUjets");
produces<reco::JetTagCollection>();
}


Expand All @@ -119,19 +123,19 @@ PixelJetPuId::~PixelJetPuId() {}

void
PixelJetPuId::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag> ("jets",edm::InputTag("hltCaloJetL1FastJetCorrected"));
desc.add<edm::InputTag> ("tracks",edm::InputTag("hltPixelTracksNoPU"));
desc.add<edm::InputTag> ("primaryVertex",edm::InputTag("hltFastPVPixelVertices"));
desc.add<double>("MinGoodJetTrackPtRatio",0.045);
desc.add<double>("MinGoodJetTrackPt",1.8);
desc.add<double>("MaxTrackDistanceToJet",0.04);
desc.add<double>("MinTrackPt",0.6);
desc.add<double>("MaxTrackChi2",20.);
desc.add<bool>("UseForwardJetsAsNoPU",true);
desc.add<double>("MinEtaForwardJets",2.4);
desc.add<double>("MinEtForwardJets",40.);
descriptions.add("pixelJetPuId",desc);
edm::ParameterSetDescription desc;
desc.add<edm::InputTag> ("jets",edm::InputTag("hltCaloJetL1FastJetCorrected"));
desc.add<edm::InputTag> ("tracks",edm::InputTag("hltPixelTracksNoPU"));
desc.add<edm::InputTag> ("primaryVertex",edm::InputTag("hltFastPVPixelVertices"));
desc.add<double>("MinGoodJetTrackPtRatio",0.045);
desc.add<double>("MinGoodJetTrackPt",1.8);
desc.add<double>("MaxTrackDistanceToJet",0.04);
desc.add<double>("MinTrackPt",0.6);
desc.add<double>("MaxTrackChi2",20.);
desc.add<bool>("UseForwardJetsAsNoPU",true);
desc.add<double>("MinEtaForwardJets",2.4);
desc.add<double>("MinEtForwardJets",40.);
descriptions.add("pixelJetPuId",desc);
}

//
Expand All @@ -141,81 +145,95 @@ PixelJetPuId::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
// ------------ method called on each new Event ------------
void PixelJetPuId::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const
{
using namespace edm;
std::unique_ptr<std::vector<reco::CaloJet> > pOut(new std::vector<reco::CaloJet> );
std::unique_ptr<std::vector<reco::CaloJet> > pOut_PUjets(new std::vector<reco::CaloJet> );

//get tracks
Handle<std::vector<reco::Track> > tracks;
iEvent.getByToken(tracksToken, tracks);
unsigned int tsize = tracks->size();
float teta[tsize], tphi[tsize];
unsigned int i=0;
for (auto const & tr : *tracks) { teta[i]=tr.eta(); tphi[i]=tr.phi();++i;}

//get jets
Handle<edm::View<reco::CaloJet> > jets;
iEvent.getByToken(jetsToken, jets);

//get primary vertices
Handle<reco::VertexCollection> primaryVertex;
iEvent.getByToken(primaryVertexToken, primaryVertex);

//get Transient Track Builder
edm::ESHandle<TransientTrackBuilder> builder;
iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);

//loop on trackIPTagInfos
if(primaryVertex->size()>0)
using namespace edm;
std::unique_ptr<std::vector<reco::CaloJet> > pOut(new std::vector<reco::CaloJet> );
std::unique_ptr<std::vector<reco::CaloJet> > pOut_PUjets(new std::vector<reco::CaloJet> );
std::unique_ptr<reco::JetTagCollection > pOut_jetTagCollection(new reco::JetTagCollection );

//get tracks
Handle<std::vector<reco::Track> > tracks;
iEvent.getByToken(tracksToken, tracks);
unsigned int tsize = tracks->size();
float teta[tsize], tphi[tsize];
unsigned int i=0;
for (auto const & tr : *tracks) { teta[i]=tr.eta(); tphi[i]=tr.phi();++i;}

//get jets
Handle<edm::View<reco::CaloJet> > jets;
iEvent.getByToken(jetsToken, jets);

Handle<edm::View<reco::Jet> > generaljets;
iEvent.getByToken(generaljetsToken, generaljets);

//get primary vertices
Handle<reco::VertexCollection> primaryVertex;
iEvent.getByToken(primaryVertexToken, primaryVertex);

//get Transient Track Builder
edm::ESHandle<TransientTrackBuilder> builder;
iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);

//init JetTagCollection
if(generaljets.product()->size()>0)
{
const reco::Vertex* pv = &*primaryVertex->begin();
//loop on jets
for(edm::View<reco::CaloJet>::const_iterator itJet = jets->begin(); itJet != jets->end(); itJet++ ) {

math::XYZVector jetMomentum = itJet->momentum();
GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());

math::XYZVector trMomentum;

//loop on tracks
if(fabs(itJet->eta())>m_mineta_fwjets)
{
if((m_fwjets) && (itJet->et()>m_minet_fwjets))
pOut->push_back(*itJet);// fill forward jet as signal jet
}
else
{
std::vector<reco::Track>::const_iterator itTrack = tracks->begin();
for (unsigned int i=0; i<tsize; ++i) {
float deltaR2=reco::deltaR2(itJet->eta(),itJet->phi(), teta[i],tphi[i]);
if(deltaR2<0.25) {
reco::TransientTrack transientTrack = builder->build(*itTrack);
float jetTrackDistance = -((IPTools::jetTrackDistance(transientTrack, direction, *pv)).second).value();

//select the tracks compabible with the jet
if(( itTrack->pt() > m_MinTrackPt) && ( itTrack->normalizedChi2() < m_MaxTrackChi2) && (jetTrackDistance<m_MaxTrackDistanceToJet))
{
trMomentum += itTrack->momentum(); //calculate the Sum(trackPt)
}
}
itTrack++;
}
//if Sum(comp.trackPt)/CaloJetPt > minPtRatio or Sum(trackPt) > minPt the jet is a signal jet
if(trMomentum.rho()/jetMomentum.rho() > m_MinGoodJetTrackPtRatio || trMomentum.rho() > m_MinGoodJetTrackPt )
{
pOut->push_back(*itJet); // fill it as signal jet
}
else//else it is a PUjet
{
pOut_PUjets->push_back(*itJet); // fill it as PUjets
}
}
}
edm::RefToBase<reco::Jet> jj = edm::RefToBase<reco::Jet>(generaljets,0);
pOut_jetTagCollection.reset(new reco::JetTagCollection(edm::makeRefToBaseProdFrom(jj, iEvent)));
}
iEvent.put(std::move(pOut));
iEvent.put(std::move(pOut_PUjets),"PUjets");


//loop on trackIPTagInfos
if(primaryVertex->size()>0)
{
const reco::Vertex* pv = &*primaryVertex->begin();
//loop on jets
for(edm::View<reco::CaloJet>::const_iterator itJet = jets->begin(); itJet != jets->end(); itJet++ ) {

math::XYZVector jetMomentum = itJet->momentum();
GlobalVector direction(jetMomentum.x(), jetMomentum.y(), jetMomentum.z());

math::XYZVector trMomentum;

if(fabs(itJet->eta())>m_mineta_fwjets)
{
if((m_fwjets) && (itJet->et()>m_minet_fwjets))
pOut->push_back(*itJet);// fill forward jet as signal jet
}
else
{
//loop on tracks
std::vector<reco::Track>::const_iterator itTrack = tracks->begin();
for (unsigned int i=0; i<tsize; ++i) {
float deltaR2=reco::deltaR2(itJet->eta(),itJet->phi(), teta[i],tphi[i]);
if(deltaR2<0.25) {
reco::TransientTrack transientTrack = builder->build(*itTrack);
float jetTrackDistance = -((IPTools::jetTrackDistance(transientTrack, direction, *pv)).second).value();

//select the tracks compabible with the jet
if(( itTrack->pt() > m_MinTrackPt) && ( itTrack->normalizedChi2() < m_MaxTrackChi2) && (jetTrackDistance<m_MaxTrackDistanceToJet))
{
trMomentum += itTrack->momentum(); //calculate the Sum(trackPt)
}
}
itTrack++;
}
//if Sum(comp.trackPt)/CaloJetPt > minPtRatio or Sum(trackPt) > minPt the jet is a signal jet
if(trMomentum.rho()/jetMomentum.rho() > m_MinGoodJetTrackPtRatio || trMomentum.rho() > m_MinGoodJetTrackPt )
{
pOut->push_back(*itJet); // fill it as signal jet
}
else//else it is a PUjet
{
pOut_PUjets->push_back(*itJet); // fill it as PUjets
}
}
RefToBase<reco::Jet> jRef(generaljets, itJet-jets->begin());
(*pOut_jetTagCollection)[jRef] = trMomentum.rho(); // fill jetTagCollection
}
}
iEvent.put(std::move(pOut));
iEvent.put(std::move(pOut_PUjets),"PUjets");
iEvent.put(std::move(pOut_jetTagCollection));
}

//define this as a plug-in
DEFINE_FWK_MODULE(PixelJetPuId);

0 comments on commit 5fb74b0

Please sign in to comment.