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 1 commit
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
230 changes: 230 additions & 0 deletions HLTrigger/Muon/interface/HLTTriMuonIsolation.h
View file
Open in desktop
@@ -0,0 +1,230 @@
#ifndef HLTrigger_Muon_HLTTriMuonIsolation_h
#define HLTrigger_Muon_HLTTriMuonIsolation_h

#include <iostream>
#include <string>

#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){
// Create the 3-muon candidates
// loop over L3/Trk muons and create all combinations
for (unsigned int i = 0; i != AllMuCands->size(); ++i){
const reco::TrackRef &tk_i = (*AllMuCands)[i].track();
Copy link
Contributor

Choose a reason for hiding this comment

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

Use some modern-style loop constructs such as for (auto const & muon : AllMuCands) {...
for these three loops.

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 think I prefer this style, as it is more legible.
auto obfuscates things, IMO, and I'm not sure that nested range-based loops would improve the legibility.

But let me try anyways and see if it ends up nicely enough.

Copy link
Contributor

Choose a reason for hiding this comment

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

I don't think you can use a range-based for loop here, since you use the index of the outer loop as the starting point of the inner loop, and so on.

It might be slightly more efficient using iterators, but the difference is probably negligible:

auto AllMuCands_end = AllMuCands->end();
for (auto i = AllMuCands->begin(); i != AllMuCands_end; ++i) {
    const reco::TrackRef &tk_i = i->track();
    for (auto j = i+1; j != AllMuCands_end; ++j) {
        const reco::TrackRef &tk_j = j->track();
        ...

and so on

Copy link
Contributor

Choose a reason for hiding this comment

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

ah yes, of course!

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, understood.
So I give up my attempts at playing magic with ranges and move to iterators

Copy link
Contributor

Choose a reason for hiding this comment

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

Or just leave the indices here, if you find it more readable. The difference is just a multiplication and a sum.

for (unsigned int j = i+1; j != AllMuCands->size(); ++j){
const reco::TrackRef &tk_j = (*AllMuCands)[j].track();
for (unsigned int k = j+1; k != AllMuCands->size(); ++k){
const reco::TrackRef &tk_k = (*AllMuCands)[k].track();

int passingPreviousFilter_1 = 0;
Copy link
Contributor

Choose a reason for hiding this comment

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

Could be calculated in the first loop

int passingPreviousFilter_2 = 0;
Copy link
Contributor

Choose a reason for hiding this comment

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

Could be calculated and kept in the 2nd loop

int passingPreviousFilter_3 = 0;
Copy link
Contributor

Choose a reason for hiding this comment

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

Only this would need to be calculated within the 3rd loop


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_i->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_1++;
if (reco::deltaR2(tk_j->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_2++;
if (reco::deltaR2(tk_k->momentum(), candTrkRef->momentum()) < (0.03*0.03)) passingPreviousFilter_3++;
}
// at least two muons must have passed the previous di-muon filter
if (((passingPreviousFilter_1 < 1) & (passingPreviousFilter_2 < 1)) ||
((passingPreviousFilter_1 < 1) & (passingPreviousFilter_3 < 1)) ||
((passingPreviousFilter_2 < 1) & (passingPreviousFilter_3 < 1))) 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((*AllMuCands)[i]);
Daughters.push_back((*AllMuCands)[j]);
Daughters.push_back((*AllMuCands)[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

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 (abs(itau->eta()) > TriMuonEtaCut_ ) continue;
Copy link
Contributor

Choose a reason for hiding this comment

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

Use std::abs instead of just abs everywhere

Copy link
Contributor

Choose a reason for hiding this comment

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

indeed, abs will truncate to an integer !

Copy link
Contributor Author

Choose a reason for hiding this comment

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

yeah, pretty lame of me

if (itau->daughter(0)->pt() < Muon1PtCut_) continue;
if (itau->daughter(1)->pt() < Muon2PtCut_) continue;
if (itau->daughter(2)->pt() < Muon3PtCut_) continue;
if ((abs(itau->charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0)) continue;
if (abs(itau->daughter(0)->vz() - itau->vz()) > MaxDZ_) continue;
if (abs(itau->daughter(1)->vz() - itau->vz()) > MaxDZ_) continue;
if (abs(itau->daughter(2)->vz() - itau->vz()) > MaxDZ_) continue;

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

// 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 (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);