btau/plugins/HLTmumutkFilter.cc

0001 #include <algorithm>
0002 #include <cmath>
0003 #include <vector>
0004
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007
0008 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0009 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0010 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0011 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0012 #include "DataFormats/VertexReco/interface/Vertex.h"
0013 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0014
0015 #include "HLTmumutkFilter.h"
0016 #include "TMath.h"
0017
0018 using namespace edm;
0019 using namespace reco;
0020 using namespace std;
0021 using namespace trigger;
0022
0023 // ----------------------------------------------------------------------
0024 HLTmumutkFilter::HLTmumutkFilter(const edm::ParameterSet& iConfig)
0025     : HLTFilter(iConfig),
0026       muCandTag_(iConfig.getParameter<edm::InputTag>("MuonTag")),
0027       muCandToken_(consumes<reco::RecoChargedCandidateCollection>(muCandTag_)),
0028       trkCandTag_(iConfig.getParameter<edm::InputTag>("TrackTag")),
0029       trkCandToken_(consumes<reco::RecoChargedCandidateCollection>(trkCandTag_)),
0030       MuMuTkVertexTag_(iConfig.getParameter<edm::InputTag>("MuMuTkVertexTag")),
0031       MuMuTkVertexToken_(consumes<reco::VertexCollection>(MuMuTkVertexTag_)),
0032       beamSpotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0033       beamSpotToken_(consumes<reco::BeamSpot>(beamSpotTag_)),
0034       maxEta_(iConfig.getParameter<double>("MaxEta")),
0035       minPt_(iConfig.getParameter<double>("MinPt")),
0036       maxNormalisedChi2_(iConfig.getParameter<double>("MaxNormalisedChi2")),
0037       minVtxProbability_(iConfig.getParameter<double>("MinVtxProbability")),
0038       minLxySignificance_(iConfig.getParameter<double>("MinLxySignificance")),
0039       minCosinePointingAngle_(iConfig.getParameter<double>("MinCosinePointingAngle")) {}
0040
0041 // ----------------------------------------------------------------------
0042 HLTmumutkFilter::~HLTmumutkFilter() = default;
0043
0044 void HLTmumutkFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0045   edm::ParameterSetDescription desc;
0046   makeHLTFilterDescription(desc);
0047   desc.add<double>("MaxEta", 2.5);
0048   desc.add<double>("MinPt", 0.0);
0049   desc.add<double>("MaxNormalisedChi2", 10.0);
0050   desc.add<double>("MinVtxProbability", 0.0);
0051   desc.add<double>("MinLxySignificance", 3.0);
0052   desc.add<double>("MinCosinePointingAngle", 0.9);
0053   desc.add<edm::InputTag>("MuonTag", edm::InputTag("hltL3MuonCandidates"));
0054   desc.add<edm::InputTag>("TrackTag", edm::InputTag("hltMumukAllConeTracks"));
0055   desc.add<edm::InputTag>("MuMuTkVertexTag", edm::InputTag("hltDisplacedmumuVtxProducerDoubleMu4Jpsi"));
0056   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOnineBeamSpot"));
0057   descriptions.add("HLTmumutkFilter", desc);
0058 }
0059
0060 // ----------------------------------------------------------------------
0061 bool HLTmumutkFilter::hltFilter(edm::Event& iEvent,
0062                                 const edm::EventSetup& iSetup,
0063                                 trigger::TriggerFilterObjectWithRefs& filterproduct) const {
0064   //get the beamspot position
0065   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0066   iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
0067   const reco::BeamSpot& vertexBeamSpot = *recoBeamSpotHandle;
0068
0069   // get vertices
0070   reco::VertexCollection displacedVertexColl;
0071   edm::Handle<reco::VertexCollection> displacedVertexCollHandle;
0072   bool foundVertexColl = iEvent.getByToken(MuMuTkVertexToken_, displacedVertexCollHandle);
0073   if (foundVertexColl)
0074     displacedVertexColl = *displacedVertexCollHandle;
0075
0076   // get muon collection
0077   Handle<RecoChargedCandidateCollection> mucands;
0078   iEvent.getByToken(muCandToken_, mucands);
0079
0080   // get track candidates around displaced muons
0081   Handle<RecoChargedCandidateCollection> trkcands;
0082   iEvent.getByToken(trkCandToken_, trkcands);
0083
0084   // Ref to Candidate object to be recorded in filter object
0085   RecoChargedCandidateRef refMu1;
0086   RecoChargedCandidateRef refMu2;
0087   RecoChargedCandidateRef refTrk;
0088
0089   if (saveTags()) {
0090     filterproduct.addCollectionTag(muCandTag_);
0091     filterproduct.addCollectionTag(trkCandTag_);
0092   }
0093
0094   bool triggered = false;
0095
0096   // loop over vertex collection
0097   reco::VertexCollection::iterator it;
0098   for (it = displacedVertexColl.begin(); it != displacedVertexColl.end(); it++) {
0099     reco::Vertex displacedVertex = *it;
0100
0101     // check if the vertex actually consists of exactly two muon + 1 track, throw exception if not
0102     if (displacedVertex.tracksSize() != 3)
0103       throw cms::Exception("BadLogic") << "HLTmumutkFilter: ERROR: the Jpsi+trk vertex must have "
0104                                        << "exactly two muons + 1 trk by definition. It now has n trakcs = "
0105                                        << displacedVertex.tracksSize() << std::endl;
0106
0107     float normChi2 = displacedVertex.normalizedChi2();
0108     if (normChi2 > maxNormalisedChi2_)
0109       continue;
0110
0111     double vtxProb = 0.0;
0112     if ((displacedVertex.chi2() >= 0.0) && (displacedVertex.ndof() > 0))
0113       vtxProb = TMath::Prob(displacedVertex.chi2(), displacedVertex.ndof());
0114     if (vtxProb < minVtxProbability_)
0115       continue;
0116
0117     // get the three tracks from the vertex
0118     auto trackIt = displacedVertex.tracks_begin();
0119     reco::TrackRef vertextkRef1 = (*trackIt).castTo<reco::TrackRef>();
0120     trackIt++;
0121     reco::TrackRef vertextkRef2 = (*trackIt).castTo<reco::TrackRef>();
0122     trackIt++;
0123     reco::TrackRef vertextkRef3 = (*trackIt).castTo<reco::TrackRef>();
0124
0125     // first find the two muon tracks in the muon collection
0126     reco::RecoChargedCandidateCollection::const_iterator mucand1;
0127     reco::RecoChargedCandidateCollection::const_iterator mucand2;
0128     reco::RecoChargedCandidateCollection::const_iterator tkcand;
0129
0130     int iFoundRefs = 0;
0131     bool track1Matched = false;
0132     bool track2Matched = false;
0133     bool track3Matched = false;
0134     for (auto cand = mucands->begin(); cand != mucands->end(); cand++) {
0135       reco::TrackRef tkRef = cand->get<reco::TrackRef>();
0136       if (!track1Matched) {
0137         if (tkRef == vertextkRef1 && iFoundRefs == 0) {
0138           mucand1 = cand;
0139           iFoundRefs++;
0140           track1Matched = true;
0141         } else if (tkRef == vertextkRef1 && iFoundRefs == 1) {
0142           mucand2 = cand;
0143           iFoundRefs++;
0144           track1Matched = true;
0145         }
0146       }
0147       if (!track2Matched) {
0148         if (tkRef == vertextkRef2 && iFoundRefs == 0) {
0149           mucand1 = cand;
0150           iFoundRefs++;
0151           track2Matched = true;
0152         } else if (tkRef == vertextkRef2 && iFoundRefs == 1) {
0153           mucand2 = cand;
0154           iFoundRefs++;
0155           track2Matched = true;
0156         }
0157       }
0158       if (!track3Matched) {
0159         if (tkRef == vertextkRef3 && iFoundRefs == 0) {
0160           mucand1 = cand;
0161           iFoundRefs++;
0162           track3Matched = true;
0163         } else if (tkRef == vertextkRef3 && iFoundRefs == 1) {
0164           mucand2 = cand;
0165           iFoundRefs++;
0166           track3Matched = true;
0167         }
0168       }
0169     }
0170     if (iFoundRefs < 2)
0171       throw cms::Exception("BadLogic") << "HLTmumutkFilterr: ERROR: the vertex must have "
0172                                        << " at least two muons by definition." << std::endl;
0173
0174     bool notYetFoundTrackRefs(true);
0175     for (auto cand = trkcands->begin(); cand != trkcands->end(); cand++) {
0176       reco::TrackRef tkRef = cand->get<reco::TrackRef>();
0177       if ((tkRef == vertextkRef1 && !track1Matched) || (tkRef == vertextkRef2 && !track2Matched) ||
0178           (tkRef == vertextkRef3 && !track3Matched)) {
0179         notYetFoundTrackRefs = false;
0180         tkcand = cand;
0181         break;
0182       }
0183     }
0184     if (notYetFoundTrackRefs)
0185       throw cms::Exception("BadLogic") << "HLTmumutkFilterr: ERROR: the vertex must have "
0186                                        << " at least one track by definition." << std::endl;
0187
0188     // calculate three-track transverse momentum
0189     math::XYZVector pperp(
0190         mucand1->px() + mucand2->px() + tkcand->px(), mucand1->py() + mucand2->py() + tkcand->py(), 0.);
0191
0192     // get vertex position and error to calculate the decay length significance
0193     const reco::Vertex::Point& vpoint = displacedVertex.position();
0194     reco::Vertex::Error verr = displacedVertex.error();
0195     GlobalPoint secondaryVertex(vpoint.x(), vpoint.y(), vpoint.z());
0196     GlobalError err(verr.At(0, 0), verr.At(1, 0), verr.At(1, 1), verr.At(2, 0), verr.At(2, 1), verr.At(2, 2));
0197
0198     GlobalPoint displacementFromBeamspot(-1 * ((vertexBeamSpot.x0() - secondaryVertex.x()) +
0199                                                (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dxdz()),
0200                                          -1 * ((vertexBeamSpot.y0() - secondaryVertex.y()) +
0201                                                (secondaryVertex.z() - vertexBeamSpot.z0()) * vertexBeamSpot.dydz()),
0202                                          0);
0203     float lxy = displacementFromBeamspot.perp();
0204     float lxyerr = sqrt(err.rerr(displacementFromBeamspot));
0205
0206     //calculate the angle between the decay length and the mumu momentum
0207     Vertex::Point vperp(displacementFromBeamspot.x(), displacementFromBeamspot.y(), 0.);
0208     float cosAlpha = vperp.Dot(pperp) / (vperp.R() * pperp.R());
0209
0210     if (pperp.R() < minPt_)
0211       continue;
0212     if (lxy / lxyerr < minLxySignificance_)
0213       continue;
0214     if (cosAlpha < minCosinePointingAngle_)
0215       continue;
0216     triggered = true;
0217
0218     refMu1 = RecoChargedCandidateRef(Ref<RecoChargedCandidateCollection>(mucands, distance(mucands->begin(), mucand1)));
0219     filterproduct.addObject(TriggerMuon, refMu1);
0220     refMu2 = RecoChargedCandidateRef(Ref<RecoChargedCandidateCollection>(mucands, distance(mucands->begin(), mucand2)));
0221     filterproduct.addObject(TriggerMuon, refMu2);
0222     refTrk =
0223         RecoChargedCandidateRef(Ref<RecoChargedCandidateCollection>(trkcands, distance(trkcands->begin(), tkcand)));
0224     filterproduct.addObject(TriggerTrack, refTrk);
0225
0226   }  //end loop vertices
0227
0228   LogDebug("HLTDisplacedMumuTrkFilter") << " >>>>> Result of HLTDisplacedMuMuTrkFilter is " << triggered;
0229   return triggered;
0230 }
0231
0232 bool HLTmumutkFilter::triggerdByPreviousLevel(const reco::RecoChargedCandidateRef& candref,
0233                                               const std::vector<reco::RecoChargedCandidateRef>& vcands) {
0234   unsigned int i = 0;
0235   unsigned int i_max = vcands.size();
0236   for (; i != i_max; ++i) {
0237     if (candref == vcands[i])
0238       return true;
0239   }
0240   return false;
0241 }