HLTTriMuonIsolation.h

Go to the documentation of this file.

#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() override;
        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           ;

        template<typename T>
        static bool ptComparer(const T & cand_1, const T & cand_2) { return cand_1.pt() > cand_2.pt(); }
        
        const double TwiceMuonMass_ = 2. * 0.1056583715; // in GeV

        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 MatchingConeSize_;
        const double MinTriMuonMass_  ;
        const double MaxTriMuonMass_  ;
        const double MaxTriMuonRadius_;
        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"       )) ,
    MatchingConeSize_                                                    (iConfig.getParameter<double>       ("MatchingConeSize"  )) ,
    MinTriMuonMass_                                                      (iConfig.getParameter<double>       ("MinTriMuonMass"    )) ,
    MaxTriMuonMass_                                                      (iConfig.getParameter<double>       ("MaxTriMuonMass"    )) , 
    MaxTriMuonRadius_                                                    (iConfig.getParameter<double>       ("MaxTriMuonRadius"  )) , 
    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-2; ++i) {
            // check that muon_i passes the previous filter
            bool passingPreviousFilter_1 = false;
            for (const auto & imu : PassedL3Muons){
                if (reco::deltaR2(i->momentum(), imu->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_1 = true;
            }
            for (auto j = i+1; j != AllMuCands_end-1; ++j) {
                // check that muon_j passes the previous filter
                bool passingPreviousFilter_2 = false;
                for (const auto & jmu : PassedL3Muons){
                    if (reco::deltaR2(j->momentum(), jmu->momentum()) < (MatchingConeSize_*MatchingConeSize_)) passingPreviousFilter_2 = true;
                }
                // if, at this point, no muons passed the previous filter just skip to the next iteration
                if (!(passingPreviousFilter_1 || passingPreviousFilter_2)) continue;
                for (auto k = j+1; k != AllMuCands_end; ++k){
                    // check that muon_k passes the previous filter
                    bool passingPreviousFilter_3 = false;
                    for (const auto & kmu : PassedL3Muons){
                        if (reco::deltaR2(k->momentum(), kmu->momentum()) < (MatchingConeSize_*MatchingConeSize_)) 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;

                    // 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;
                    daughters.reserve(3);
                    
                    daughters.push_back(*i);
                    daughters.push_back(*j);
                    daughters.push_back(*k);
                                                            
                    std::sort(daughters.begin(), daughters.end(), ptComparer<reco::RecoChargedCandidate>);

                    // Muon kinematic selections
                    if (daughters[0].pt() < Muon1PtCut_) continue;
                    if (daughters[1].pt() < Muon2PtCut_) continue;
                    if (daughters[2].pt() < Muon3PtCut_) continue;

                    // assign the tau its daughters
                    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 to each other 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;
                    
                    // require muons to be collimated
                    bool collimated = true;           
                    for (auto const &idau : daughters){
                        if (reco::deltaR2(tau.p4(), idau.p4()) > MaxTriMuonRadius_*MaxTriMuonRadius_) {
                            collimated = false;
                            break;
                        }
                    }
                                  
                    if (!collimated) continue;
                                            
                    // Tau kinematic selections
                    if (tau.pt()   < TriMuonPtCut_  ) continue;
                    if (tau.mass() < MinTriMuonMass_) continue;
                    if (tau.mass() > MaxTriMuonMass_) continue;
                    if (std::abs(tau.eta()) > TriMuonEtaCut_ ) continue;

                    // Tau charge selection 
                    if ((std::abs(tau.charge()) != TriMuonAbsCharge_) & (TriMuonAbsCharge_ >= 0)) continue;
 
                    // Sanity check against duplicates, di-muon masses must be > 2 * mass_mu
                    if ( (tau.daughter(0)->p4() + tau.daughter(1)->p4()).mass() < TwiceMuonMass_) continue;
                    if ( (tau.daughter(0)->p4() + tau.daughter(2)->p4()).mass() < TwiceMuonMass_) continue;
                    if ( (tau.daughter(1)->p4() + tau.daughter(2)->p4()).mass() < TwiceMuonMass_) continue;
                    
                    // a good tau, at last                      
                    Taus->push_back(tau);
                }
            }
        }

        // Sort taus by pt
        std::sort(Taus->begin(), Taus->end(), ptComparer<reco::CompositeCandidate>);

        // Loop over taus and further select by isolation
        for (const auto & itau : *Taus){

            // remove the candidate pt from the iso sum
            double sumPt = -itau.pt();
            
            // compute iso sum pT
            for (const auto & itrk : *IsoTracks){
                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
            double chAbsIsoCut = EnableAbsIso_ ? ChargedAbsIsoCut_             : std::numeric_limits<double>::infinity();
            double chRelIsoCut = EnableRelIso_ ? ChargedRelIsoCut_ * itau.pt() : std::numeric_limits<double>::infinity();
            
            if ( !( (sumPt < chAbsIsoCut)||(sumPt < chRelIsoCut) ) ) 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>("MatchingConeSize", 0.03 );
  desc.add<double>("MinTriMuonMass"  , 0.5  );
  desc.add<double>("MaxTriMuonMass"  , 2.8  );
  desc.add<double>("MaxTriMuonRadius", 0.6  );
  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