Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

L1 pf jets matching - Adding the template module #19068

Merged
merged 9 commits into from Jun 8, 2017
207 changes: 207 additions & 0 deletions RecoTauTag/HLTProducers/interface/L1TJetsMatching.h
@@ -0,0 +1,207 @@
// -*- C++ -*-
//
// Package: RecoTauTag/HLTProducers
// Class: L1TJetsMatching
//
/**\class L1TJetsMatching L1TJetsMatching.h
RecoTauTag/HLTProducers/interface/L1TJetsMatching.h
Description:
Matching L1 to PF/Calo Jets. Used for HLT_VBF paths.
*Matches PF/Calo Jets to L1 jets from the dedicated seed
*Adds selection criteria to the leading/subleading jets as well as the maximum dijet mass
*Separates collections of PF/Calo jets into two categories


*/
//
// Original Author: Vukasin Milosevic
// Created: Thu, 01 Jun 2017 17:23:00 GMT
//
//




#ifndef RecoTauTag_HLTProducers_L1TJetsMatching_h
#define RecoTauTag_HLTProducers_L1TJetsMatching_h

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/JetReco/interface/PFJetCollection.h"
#include "DataFormats/TauReco/interface/PFTauFwd.h"
#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"

#include "Math/GenVector/VectorUtil.h"
#include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "DataFormats/JetReco/interface/PFJet.h"

#include "HLTrigger/HLTcore/interface/defaultModuleLabel.h"
#include "DataFormats/Math/interface/deltaR.h"

#include <map>
#include <vector>

template< typename T>
class L1TJetsMatching: public edm::global::EDProducer<> {
public:
explicit L1TJetsMatching(const edm::ParameterSet&);
~L1TJetsMatching();
virtual void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:

const edm::EDGetTokenT<std::vector<T>> jetSrc_;
const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> jetTrigger_;
const double pt1_Min;
const double pt2_Min;
const double mjj_Min;

};
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Identify all class data members by a trailing underscore "_" such as pt1Min_ etc.

//
// class decleration
//
using namespace reco ;
using namespace std ;
using namespace edm ;
using namespace trigger;


template< typename T>
std::pair<std::vector<T>,std::vector<T>> categorise(const std::vector<T>& pfMatchedJets, double pt1, double pt2, double Mjj)
{
std::pair<std::vector<T>,std::vector<T>> output;
unsigned int i1 = 0;
unsigned int i2 = 0;
double mjj = 0;
if (pfMatchedJets.size()>1){
for (unsigned int i = 0; i < pfMatchedJets.size()-1; i++){

const T & myJet1 = (pfMatchedJets)[i];

for (unsigned int j = i+1; j < pfMatchedJets.size(); j++)
{
const T & myJet2 = (pfMatchedJets)[j];

const double mjj_test = (myJet1.p4()+myJet2.p4()).M();

if (mjj_test > mjj){

mjj =mjj_test;
i1 = i;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i1 and i2 are not used. they can be removed.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i1 and i2 are used in remembering the indices of the two jets forming the largest mjj

i2 = j;
}
}
}

const T & myJet1 = (pfMatchedJets)[i1];
const T & myJet2 = (pfMatchedJets)[i2];

if ((mjj > Mjj) && (myJet1.pt() >= pt1) && (myJet2.pt() > pt2) )
{

output.first.push_back(myJet1);
output.first.push_back(myJet2);

}

if ((mjj > Mjj) && (myJet1.pt() < pt1) && (myJet1.pt() > pt2) && (myJet2.pt() > pt2))
{
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

myJet1.pt()>pt2? Why pt2?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the three jets case, where the two jets forming the max mjj have to have pt more than 35GeV but the higher pt one (myJets1) has pt<95 (otherwise it would have been the above two jets case).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I should know but I do not - the jets in your pfMatchedJets vector are increasing pt ordered?


const T & myJetTest = (pfMatchedJets)[0];
if (myJetTest.pt()>pt1){
output.second.push_back(myJet1);
output.second.push_back(myJet2);
output.second.push_back(myJetTest);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do you need to protect against i1==0 || i2==0?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i2 should not be able to be ==0 by construction (j=i+1,...). Also, I think that it is safe even without i1==0 protection, because of the pt>pt1 requirement (in the if above there is a req for myJet1.pt<=pt1 so if i1 = 0 it will fail the statement at the end)


}
}

}

return output;

}
template< typename T>
L1TJetsMatching<T>::L1TJetsMatching(const edm::ParameterSet& iConfig):
jetSrc_ ( consumes<std::vector<T>> (iConfig.getParameter<InputTag>("JetSrc" ) ) ),
jetTrigger_( consumes<trigger::TriggerFilterObjectWithRefs>(iConfig.getParameter<InputTag>("L1JetTrigger") ) ),
pt1_Min ( iConfig.getParameter<double>("pt1_Min")),
pt2_Min ( iConfig.getParameter<double>("pt2_Min")),
mjj_Min ( iConfig.getParameter<double>("mjj_Min"))
{
produces<std::vector<T>>("TwoJets");
produces<std::vector<T>>("ThreeJets");
}
template< typename T>
L1TJetsMatching<T>::~L1TJetsMatching(){ }

template< typename T>
void L1TJetsMatching<T>::produce(edm::StreamID iSId, edm::Event& iEvent, const edm::EventSetup& iES) const
{

unique_ptr<std::vector<T>> pfMatchedJets(new std::vector<T>);
std::pair<std::vector<T>,std::vector<T>> output;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could you not avoid the copy of the output jets by making your interface a

std::pair< unique_ptr<std::vector>, unique_ptr<std::vector> >

?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure if I understand the question.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you currently do

std::pair<vector,vector> output = categorize(...);
unique_ptr<std::vector> output1(new std::vector(output.first));
unique_ptr<std::vector> output2(new std::vector(output.second));
iEvent.put(std::move(output1),"TwoJets");
iEvent.put(std::move(output2),"ThreeJets");

[eg, you create each vector 3 times]

instead something like

std::pair< unique_ptr<std::vector>, unique_ptr<std::vector> > output = categorize(...);
iEvent.put(std::move(output.first),"TwoJets");
iEvent.put(std::move(output.second),"ThreeJets");

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried doing similar thing at the very beginning but couldn't make it work with doing .put(output.first). I can try again if it is really an issue?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

copying things is an issue - what problems did you run into?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I couldn't make it compile at the beginning so my colleagues showed me how they were doing it in their analysis code. Do I need to change this?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lets fix as a followup PR this week. please let me know what problems you had.


double deltaR = 1.0;
double matchingR = 0.5;

// Getting HLT jets to be matched
edm::Handle<std::vector<T> > pfJets;
iEvent.getByToken( jetSrc_, pfJets );

edm::Handle<trigger::TriggerFilterObjectWithRefs> l1TriggeredJets;
iEvent.getByToken(jetTrigger_,l1TriggeredJets);

//l1t::TauVectorRef jetCandRefVec;
l1t::JetVectorRef jetCandRefVec;
l1TriggeredJets->getObjects( trigger::TriggerL1Jet,jetCandRefVec);

math::XYZPoint a(0.,0.,0.);

//std::cout<<"PFsize= "<<pfJets->size()<<endl<<" L1size= "<<jetCandRefVec.size()<<std::endl;
for(unsigned int iJet = 0; iJet < pfJets->size(); iJet++){
const T & myJet = (*pfJets)[iJet];
for(unsigned int iL1Jet = 0; iL1Jet < jetCandRefVec.size(); iL1Jet++){
// Find the relative L2pfJets, to see if it has been reconstructed
// if ((iJet<3) && (iL1Jet==0)) std::cout<<myJet.p4().Pt()<<" ";
deltaR = reco::deltaR2(myJet.p4().Vect(), (jetCandRefVec[iL1Jet]->p4()).Vect());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think deltaR2 returns deltaR^2... please check and then square the cut value (once at the beginning) and compare, rather than taking the sqrt, to be faster.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, deltaR2 is deltaR^2

if(deltaR < matchingR ) {
pfMatchedJets->push_back(myJet);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could eliminate deltaR2 by

if (reco::deltaR2(myJet.p4().Vect(), (jetCandRefVec[iL1Jet]->p4()).Vect()) < matchingR2) {

Also, I do not think you need the .vect() but could pass .p4() directly

break;
}
}
}

output= categorise(*pfMatchedJets,pt1_Min,pt2_Min, mjj_Min);
unique_ptr<std::vector<T>> output1(new std::vector<T>(output.first));
unique_ptr<std::vector<T>> output2(new std::vector<T>(output.second));

iEvent.put(std::move(output1),"TwoJets");
iEvent.put(std::move(output2),"ThreeJets");

}
template< typename T>
void L1TJetsMatching<T>::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
{
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("L1JetTrigger", edm::InputTag("hltL1DiJetVBF"))->setComment("Name of trigger filter" );
desc.add<edm::InputTag>("JetSrc" , edm::InputTag("hltAK4PFJetsTightIDCorrected"))->setComment("Input collection of PFJets");
desc.add<double> ("pt1_Min",95.0)->setComment("Minimal pT1 of PFJets to match");
desc.add<double> ("pt2_Min",35.0)->setComment("Minimal pT2 of PFJets to match");
desc.add<double> ("mjj_Min",650.0)->setComment("Minimal mjj of matched PFjets");
descriptions.setComment("This module produces collection of PFJetss matched to L1 Taus / Jets passing a HLT filter (Only p4 and vertex of returned PFJetss are set).");
descriptions.add(defaultModuleLabel<L1TJetsMatching<T>>(), desc);
}



#endif
9 changes: 9 additions & 0 deletions RecoTauTag/HLTProducers/src/SealModule.cc
Expand Up @@ -20,6 +20,13 @@
#include "HLTPFTauPairLeadTrackDzMatchFilter.h"
#include "RecoTauTag/HLTProducers/interface/L2TauPixelIsoTagProducer.h"

#include "DataFormats/JetReco/interface/PFJet.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "RecoTauTag/HLTProducers/interface/L1TJetsMatching.h"

typedef L1TJetsMatching<reco::PFJet> L1TPFJetsMatching ;
typedef L1TJetsMatching<reco::CaloJet> L1TCaloJetsMatching ;

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

L1T.... instead of L1.... (new L1T convention)

DEFINE_EDM_PLUGIN(TrackingRegionProducerFactory, TauRegionalPixelSeedGenerator, "TauRegionalPixelSeedGenerator");
DEFINE_EDM_PLUGIN(TrackingRegionProducerFactory, TrackingRegionsFromBeamSpotAndL2Tau, "TrackingRegionsFromBeamSpotAndL2Tau");
DEFINE_EDM_PLUGIN(TrackingRegionProducerFactory, CandidateSeededTrackingRegionsProducer, "CandidateSeededTrackingRegionsProducer");
Expand All @@ -45,4 +52,6 @@ DEFINE_FWK_MODULE(VertexFromTrackProducer);
//DEFINE_FWK_MODULE(L2TauPixelTrackMatch);
DEFINE_FWK_MODULE(HLTPFTauPairLeadTrackDzMatchFilter);
DEFINE_FWK_MODULE(L2TauPixelIsoTagProducer);
DEFINE_FWK_MODULE(L1TCaloJetsMatching);
DEFINE_FWK_MODULE(L1TPFJetsMatching);

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

L1T... instead of L1...