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.

Sign up for GitHub

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

Jump to bottom

Add tau3mu HLT producer #19696

Merged
merged 7 commits into from Jul 18, 2017
Merged

Add tau3mu HLT producer #19696

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e0b3aa8
add tau3mu HLT producer
rmanzoni Jul 12, 2017
3a1a02a
improve loop logic, anticipate dZ cut
rmanzoni Jul 14, 2017
049a2d3
spare a loop if you can
rmanzoni Jul 14, 2017
32439c5
remove unneeded include from boost lib
rmanzoni Jul 14, 2017
a3a2fd1
add option to require collimated muons
rmanzoni Jul 14, 2017
97140bd
pass by reference, improve loops, stick to CMS' variable naming style
rmanzoni Jul 17, 2017
dc4fe75
shorten loops when muons run out
rmanzoni Jul 17, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions HLTrigger/Muon/BuildFile.xml
View file
Open in desktop
@@ -1,3 +1,4 @@
<use name="DataFormats/Candidate"/>
<use name="DataFormats/Common"/>
<use name="DataFormats/HLTReco"/>
<use name="DataFormats/L1GlobalMuonTrigger"/>
Expand Down
242 changes: 242 additions & 0 deletions HLTrigger/Muon/interface/HLTTriMuonIsolation.h
View file
Open in desktop
@@ -0,0 +1,242 @@
#ifndef HLTrigger_Muon_HLTTriMuonIsolation_h
#define HLTrigger_Muon_HLTTriMuonIsolation_h

#include <iostream>
#include <string>
#include <boost/range/irange.hpp>

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Candidate/interface/CompositeCandidate.h"
#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "DataFormats/MuonReco/interface/MuonFwd.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
#include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"


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

private:
const edm::EDGetTokenT<reco::RecoChargedCandidateCollection> L3MuonsToken_ ;
const edm::EDGetTokenT<reco::RecoChargedCandidateCollection> AllMuonsToken_ ;
const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> L3DiMuonsFilterToken_;
const edm::EDGetTokenT<reco::TrackCollection> IsoTracksToken_ ;

edm::Handle<reco::RecoChargedCandidateCollection> L3MuCands ;
edm::Handle<trigger::TriggerFilterObjectWithRefs> L3DiMuonsFilterCands;
edm::Handle<reco::RecoChargedCandidateRef> PassedL3Muons ;
edm::Handle<reco::RecoChargedCandidateCollection> AllMuCands ;
edm::Handle<reco::TrackCollection> IsoTracks ;

static bool ptComparer(const reco::RecoChargedCandidate mu_1, const reco::RecoChargedCandidate mu_2) { return mu_1.pt() > mu_2.pt(); }
Copy link
Contributor

Choose a reason for hiding this comment

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

hi @rmanzoni - please pass by const &


const double Muon1PtCut_ ;
const double Muon2PtCut_ ;
const double Muon3PtCut_ ;
const double TriMuonPtCut_ ;
const double TriMuonEtaCut_ ;
const double ChargedRelIsoCut_;
const double ChargedAbsIsoCut_;
const double IsoConeSize_ ;
const double MinTriMuonMass_ ;
const double MaxTriMuonMass_ ;
const int TriMuonAbsCharge_;
const double MaxDZ_ ;
const bool EnableRelIso_ ;
const bool EnableAbsIso_ ;
};

HLTTriMuonIsolation::HLTTriMuonIsolation(const edm::ParameterSet& iConfig):
L3MuonsToken_ (consumes<reco::RecoChargedCandidateCollection> (iConfig.getParameter<edm::InputTag>("L3MuonsSrc" ))),
AllMuonsToken_ (consumes<reco::RecoChargedCandidateCollection> (iConfig.getParameter<edm::InputTag>("AllMuonsSrc" ))),
L3DiMuonsFilterToken_(consumes<trigger::TriggerFilterObjectWithRefs> (iConfig.getParameter<edm::InputTag>("L3DiMuonsFilterSrc"))),
IsoTracksToken_ (consumes<reco::TrackCollection> (iConfig.getParameter<edm::InputTag>("IsoTracksSrc" ))),
Muon1PtCut_ (iConfig.getParameter<double> ("Muon1PtCut" )) ,
Muon2PtCut_ (iConfig.getParameter<double> ("Muon2PtCut" )) ,
Muon3PtCut_ (iConfig.getParameter<double> ("Muon3PtCut" )) ,
TriMuonPtCut_ (iConfig.getParameter<double> ("TriMuonPtCut" )) ,
TriMuonEtaCut_ (iConfig.getParameter<double> ("TriMuonEtaCut" )) ,
ChargedRelIsoCut_ (iConfig.getParameter<double> ("ChargedRelIsoCut" )) ,
ChargedAbsIsoCut_ (iConfig.getParameter<double> ("ChargedAbsIsoCut" )) ,
IsoConeSize_ (iConfig.getParameter<double> ("IsoConeSize" )) ,
MinTriMuonMass_ (iConfig.getParameter<double> ("MinTriMuonMass" )) ,
MaxTriMuonMass_ (iConfig.getParameter<double> ("MaxTriMuonMass" )) ,
TriMuonAbsCharge_ (iConfig.getParameter<int> ("TriMuonAbsCharge" )) ,
MaxDZ_ (iConfig.getParameter<double> ("MaxDZ" )) ,
EnableRelIso_ (iConfig.getParameter<bool> ("EnableRelIso" )) ,
EnableAbsIso_ (iConfig.getParameter<bool> ("EnableAbsIso" ))
{
//register products
produces<reco::CompositeCandidateCollection>("Taus");
produces<reco::CompositeCandidateCollection>("SelectedTaus");
}

HLTTriMuonIsolation::~HLTTriMuonIsolation(){ }

void
HLTTriMuonIsolation::produce(edm::StreamID sid, edm::Event & iEvent, edm::EventSetup const & iSetup) const
{
std::unique_ptr<reco::CompositeCandidateCollection> Taus (new reco::CompositeCandidateCollection);
std::unique_ptr<reco::CompositeCandidateCollection> SelectedTaus(new reco::CompositeCandidateCollection);

// Get the L3 muon candidates
edm::Handle<reco::RecoChargedCandidateCollection> L3MuCands;
iEvent.getByToken(L3MuonsToken_, L3MuCands);

// Get the L3 muon candidates that passed the filter
edm::Handle<trigger::TriggerFilterObjectWithRefs> L3DiMuonsFilterCands;
iEvent.getByToken(L3DiMuonsFilterToken_, L3DiMuonsFilterCands);

std::vector<reco::RecoChargedCandidateRef> PassedL3Muons;
L3DiMuonsFilterCands->getObjects(trigger::TriggerMuon, PassedL3Muons);

// Get the Trk + L3 muon candidates (after merging)
edm::Handle<reco::RecoChargedCandidateCollection> AllMuCands;
iEvent.getByToken(AllMuonsToken_, AllMuCands);

// Get iso tracks
edm::Handle<reco::TrackCollection> IsoTracks;
iEvent.getByToken(IsoTracksToken_, IsoTracks);

if (AllMuCands->size() >= 3 && L3MuCands->size() >= 2){
// Create the 3-muon candidates
// loop over L3/Trk muons and create all combinations
auto AllMuCands_end = AllMuCands->end();
for (auto i = AllMuCands->begin(); i != AllMuCands_end; ++i) {
// check that muon_i passes the previous filter
bool passingPreviousFilter_1 = false;
const reco::TrackRef &tk_i = i->track();
for (std::vector<reco::RecoChargedCandidateRef>::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){
Copy link
Contributor

Choose a reason for hiding this comment

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

please make this a range-based loop

Copy link
Contributor Author

Choose a reason for hiding this comment

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

hi @davidlange6,

this point was brought up during the first review of the code.
Blame my poor skills, but I can't think of a neat way to use range based loops in this case, where I need to build all class-3 combinations and I need to know the index of the iteration in each loop.

Do you have any suggestions?

Copy link
Contributor

Choose a reason for hiding this comment

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

i didn't comment on the outer loop and the next loop because indeed they are not simple (though possible) - i had assumed these second loops can be

for ( auto const & imu : PassedL3Muons) {
reco::TrackRef candTrkRef = imu->getreco::TrackRef();
....

(but why do you need the trackref at all? btw)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

ok, you're right, I misread, sorry.

Now that it's dR based in fact I probably don't need the TkRef but I'm fine with the muon p4 alone.

The reason why it's written this way is that, in principle, the product of the filter and the muon track, would point to the same reference so one could simply check the matching by checking the reference.
The actual case is a bit more complex than that, because my muons undergo a second reconstruction and their original track might be replaced by a new one, therefore I fallback on the plain dR matching.

reco::TrackRef candTrkRef = (*imu)->get<reco::TrackRef>();
if (reco::deltaR2(tk_i->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_1 = true;
Copy link
Contributor

Choose a reason for hiding this comment

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

0.03 is a magic number - better make it a configurable plugin parameter and of course add it to fillDescriptions()

Copy link
Contributor Author

Choose a reason for hiding this comment

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

ok

}
for (auto j = i+1; j != AllMuCands_end; ++j) {
Copy link
Contributor

Choose a reason for hiding this comment

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

please protect against i=AllMuCands_end-1

// check that muon_j passes the previous filter
bool passingPreviousFilter_2 = false;
const reco::TrackRef &tk_j = j->track();
for (std::vector<reco::RecoChargedCandidateRef>::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){
Copy link
Contributor

Choose a reason for hiding this comment

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

range-based loop

reco::TrackRef candTrkRef = (*imu)->get<reco::TrackRef>();
if (reco::deltaR2(tk_j->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_2 = true;
}
for (auto k = j+1; k != AllMuCands_end; ++k){
Copy link
Contributor

Choose a reason for hiding this comment

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

No need to execute 3rd loop if passingPreviousFilter_1 and passingPreviousFilter_2 are both FALSE (or 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.

good point

// check that muon_k passes the previous filter
bool passingPreviousFilter_3 = false;
const reco::TrackRef &tk_k = k->track();
for (std::vector<reco::RecoChargedCandidateRef>::const_iterator imu = PassedL3Muons.begin(); imu != PassedL3Muons.end(); ++imu){
reco::TrackRef candTrkRef = (*imu)->get<reco::TrackRef>();
if (reco::deltaR2(tk_k->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_3 = true;
}
// at least two muons must have passed the previous di-muon filter
if (!( (passingPreviousFilter_1 & passingPreviousFilter_2 ) ||
(passingPreviousFilter_1 & passingPreviousFilter_3 ) ||
(passingPreviousFilter_2 & passingPreviousFilter_3 ) )) continue;

Copy link
Contributor

Choose a reason for hiding this comment

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

IOW, the above should be caluclated less often in the corresponding loop for i/1, j/2 and k/3.

Copy link
Contributor

Choose a reason for hiding this comment

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

Actually, if the passingPreviousFilter_1/2/3 were all unsigned ints - either 0 or 1,
then the condition would be simply passingPreviousFilter_1+passingPreviousFilter_2+passingPreviousFilter_3<3

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 it, but I found out that you can have the same muon appearing twice in the output of the previous filter, just because there is more than one combination that picks it up.

For example if the two good pairs are (1,2) and (1,3), you end up having 1 appearing twice and that would bump the sum up artificially

Copy link
Contributor Author

@rmanzoni rmanzoni Jul 14, 2017 •
edited

Choose a reason for hiding this comment

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

In light of the comment you made below, here the condition reduces to checking that
passingPreviousFilter_3 is true, correct?
wrong...

Copy link
Contributor

Choose a reason for hiding this comment

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

You mean the same muons appears twice in the collection?
Otherwise, the nested loops start at +1 index off the outer loop index, so the same
muon is not used twice.

Copy link
Contributor

Choose a reason for hiding this comment

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

IOW, even with (1,2) and (1,3), passingPreviousFilter_1 would be at most 1.

Copy link
Contributor Author

@rmanzoni rmanzoni Jul 14, 2017 •
edited

Choose a reason for hiding this comment

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

the same muon can appear more than once in PassedL3Muons, so if I loop on the entire PassedL3Muons collection for each muon of AllMuCands, even in the nested loop, that muon of AllMuCands can bump the counter several times

// Create a composite candidate to be a tau
reco::CompositeCandidate Tau;

// sort the muons by pt and add them to the tau
reco::RecoChargedCandidateCollection Daughters;
Copy link
Contributor

Choose a reason for hiding this comment

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

please use lower-case for variable names: daughters, tau


Daughters.push_back(*i);
Copy link
Contributor

Choose a reason for hiding this comment

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

reserve(3)?

Daughters.push_back(*j);
Daughters.push_back(*k);

std::sort(Daughters.begin(), Daughters.end(), ptComparer);

Tau.addDaughter((Daughters)[0], "Muon_1");
Tau.addDaughter((Daughters)[1], "Muon_2");
Tau.addDaughter((Daughters)[2], "Muon_3");

// start building the tau
int charge = Daughters[0].charge() + Daughters[1].charge() + Daughters[2].charge();
math::XYZTLorentzVectorD taup4 = Daughters[0].p4() + Daughters[1].p4() + Daughters[2].p4() ;
int tauPdgId = charge > 0? 15 : -15;

Tau.setP4(taup4);
Tau.setCharge(charge);
Tau.setPdgId(tauPdgId);
Tau.setVertex((Daughters)[0].vertex()); // assign the leading muon vertex as tau vertex

// the three muons must be close in Z
if (std::abs(Tau.daughter(0)->vz() - Tau.vz()) > MaxDZ_) continue;
if (std::abs(Tau.daughter(1)->vz() - Tau.vz()) > MaxDZ_) continue;
if (std::abs(Tau.daughter(2)->vz() - Tau.vz()) > MaxDZ_) continue;

Taus->push_back(Tau);
}
}
}

// Loop over taus and further select
for (reco::CompositeCandidateCollection::const_iterator itau = Taus->begin(); itau != Taus->end(); ++itau){
if ( itau->pt() < TriMuonPtCut_ ) continue;
if ( itau->mass() < MinTriMuonMass_) continue;
if ( itau->mass() > MaxTriMuonMass_) continue;
if (std::abs(itau->eta()) > TriMuonEtaCut_ ) continue;
if (itau->daughter(0)->pt() < Muon1PtCut_ ) continue;
if (itau->daughter(1)->pt() < Muon2PtCut_ ) continue;
if (itau->daughter(2)->pt() < Muon3PtCut_ ) continue;
if ((std::abs(itau->charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0)) continue;

// remove the candidate pt from the iso sum
double sumPt = -itau->pt();

// compute iso sum pT
for (reco::TrackCollection::const_iterator itrk = IsoTracks->begin(); itrk != IsoTracks->end(); ++itrk){
Copy link
Contributor

Choose a reason for hiding this comment

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

range-based loop will work well here

if (reco::deltaR2(itrk->momentum(), itau->p4()) > IsoConeSize_*IsoConeSize_) continue;
if (std::abs(itrk->vz() - itau->vz()) > MaxDZ_) continue;
sumPt += itrk->pt();
}

// apply the isolation cut
if ((std::max(0., sumPt) > (EnableAbsIso_ * ChargedAbsIsoCut_)) ||
Copy link
Contributor

Choose a reason for hiding this comment

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

why the std::max?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

in practice it is not needed.

I think I put it because I didn't like the idea of a negative sum pt (whilst it might happen because tracks from different tracking sequences are used at the same time).

Will remove it altogether.

(std::max(0., sumPt) > (EnableRelIso_ * ChargedRelIsoCut_ * itau->pt()))) continue;

SelectedTaus->push_back(*itau);
}
}

// finally put the vector of 3-muon candidates in the event
iEvent.put(std::move(Taus) , "Taus" );
iEvent.put(std::move(SelectedTaus), "SelectedTaus");
}

void
HLTTriMuonIsolation::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
{
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("L3MuonsSrc" , edm::InputTag("hltIterL3FromL2MuonCandidates" ));
desc.add<edm::InputTag>("AllMuonsSrc" , edm::InputTag("hltGlbTrkMuonCands" ));
desc.add<edm::InputTag>("L3DiMuonsFilterSrc", edm::InputTag("hltDiMuonForTau3MuDzFiltered0p3"));
desc.add<edm::InputTag>("IsoTracksSrc" , edm::InputTag("hltIter2L3FromL2MuonMerged" ));
desc.add<double>("Muon1PtCut" , 5. );
desc.add<double>("Muon2PtCut" , 3. );
desc.add<double>("Muon3PtCut" , 0. );
desc.add<double>("TriMuonPtCut" , 8. );
desc.add<double>("TriMuonEtaCut" , 2.5 );
desc.add<double>("ChargedAbsIsoCut", 3.0 );
desc.add<double>("ChargedRelIsoCut", 0.1 );
desc.add<double>("IsoConeSize" , 0.5 );
desc.add<double>("MinTriMuonMass" , 0.5 );
desc.add<double>("MaxTriMuonMass" , 2.8 );
desc.add<int> ("TriMuonAbsCharge", -1 );
desc.add<double>("MaxDZ" , 0.3 );
desc.add<bool> ("EnableRelIso" , false);
desc.add<bool> ("EnableAbsIso" , true );
descriptions.add("hltTriMuonIsolationProducer",desc);
}

#endif
3 changes: 3 additions & 0 deletions HLTrigger/Muon/src/SealModule.cc
View file
Open in desktop
Expand Up @@ -45,3 +45,6 @@ DEFINE_FWK_MODULE(HLTL1MuonSelector);
DEFINE_FWK_MODULE(HLTL1TMuonSelector);
DEFINE_FWK_MODULE(HLTScoutingMuonProducer);
DEFINE_FWK_MODULE(HLTL1MuonNoL2Selector);

#include "HLTrigger/Muon/interface/HLTTriMuonIsolation.h"
DEFINE_FWK_MODULE(HLTTriMuonIsolation);