btau/plugins/HLTmumutkVtxProducer.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 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0008 #include "DataFormats/Math/interface/deltaR.h"
0009 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0010 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0011 #include "DataFormats/VertexReco/interface/Vertex.h"
0012 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0013 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0014 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0015 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0016 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0017 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0018 #include "HLTmumutkVtxProducer.h"
0019
0020 using namespace edm;
0021 using namespace reco;
0022 using namespace std;
0023 using namespace trigger;
0024
0025 HLTmumutkVtxProducer::HLTmumutkVtxProducer(const edm::ParameterSet& iConfig)
0026     : transientTrackRecordToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0027       muCandTag_(iConfig.getParameter<edm::InputTag>("MuCand")),
0028       muCandToken_(consumes<reco::RecoChargedCandidateCollection>(muCandTag_)),
0029       trkCandTag_(iConfig.getParameter<edm::InputTag>("TrackCand")),
0030       trkCandToken_(consumes<reco::RecoChargedCandidateCollection>(trkCandTag_)),
0031       previousCandTag_(iConfig.getParameter<edm::InputTag>("PreviousCandTag")),
0032       previousCandToken_(consumes<trigger::TriggerFilterObjectWithRefs>(previousCandTag_)),
0033       mfName_(iConfig.getParameter<std::string>("SimpleMagneticField")),
0034       idealMagneticFieldRecordToken_(esConsumes(edm::ESInputTag("", mfName_))),
0035       thirdTrackMass_(iConfig.getParameter<double>("ThirdTrackMass")),
0036       maxEta_(iConfig.getParameter<double>("MaxEta")),
0037       minPt_(iConfig.getParameter<double>("MinPt")),
0038       minInvMass_(iConfig.getParameter<double>("MinInvMass")),
0039       maxInvMass_(iConfig.getParameter<double>("MaxInvMass")),
0040       minD0Significance_(iConfig.getParameter<double>("MinD0Significance")),
0041       // minimum delta-R^2 threshold (with sign) for non-overlapping tracks
0042       overlapDR2_(iConfig.getParameter<double>("OverlapDR") * std::abs(iConfig.getParameter<double>("OverlapDR"))),
0043       beamSpotTag_(iConfig.getParameter<edm::InputTag>("BeamSpotTag")),
0044       beamSpotToken_(consumes<reco::BeamSpot>(beamSpotTag_)) {
0045   produces<VertexCollection>();
0046 }
0047
0048 void HLTmumutkVtxProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0049   edm::ParameterSetDescription desc;
0050   desc.add<edm::InputTag>("MuCand", edm::InputTag("hltMuTracks"));
0051   desc.add<edm::InputTag>("TrackCand", edm::InputTag("hltMumukAllConeTracks"));
0052   desc.add<edm::InputTag>("PreviousCandTag", edm::InputTag("hltDisplacedmumuFilterDoubleMu4Jpsi"));
0053   desc.add<std::string>("SimpleMagneticField", "");
0054   desc.add<double>("ThirdTrackMass", 0.493677);
0055   desc.add<double>("MaxEta", 2.5);
0056   desc.add<double>("MinPt", 3.0);
0057   desc.add<double>("MinInvMass", 0.0);
0058   desc.add<double>("MaxInvMass", 99999.);
0059   desc.add<double>("MinD0Significance", 0.0);
0060   desc.add<double>("OverlapDR", 1.44e-4);
0061   desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("hltOfflineBeamSpot"));
0062   descriptions.add("HLTmumutkVtxProducer", desc);
0063 }
0064
0065 void HLTmumutkVtxProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0066   const double MuMass(0.106);
0067   const double MuMass2(MuMass * MuMass);
0068   const double thirdTrackMass2(thirdTrackMass_ * thirdTrackMass_);
0069
0070   // get hold of muon trks
0071   Handle<RecoChargedCandidateCollection> mucands;
0072   iEvent.getByToken(muCandToken_, mucands);
0073
0074   // get the transient track builder
0075   auto const& theB = iSetup.getHandle(transientTrackRecordToken_);
0076
0077   //get the beamspot position
0078   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0079   iEvent.getByToken(beamSpotToken_, recoBeamSpotHandle);
0080
0081   //get the b field
0082   auto const& bFieldHandle = iSetup.getHandle(idealMagneticFieldRecordToken_);
0083   const MagneticField* magField = bFieldHandle.product();
0084   TSCBLBuilderNoMaterial blsBuilder;
0085
0086   // get track candidates around displaced muons
0087   Handle<RecoChargedCandidateCollection> trkcands;
0088   iEvent.getByToken(trkCandToken_, trkcands);
0089
0090   unique_ptr<VertexCollection> vertexCollection(new VertexCollection());
0091
0092   // Ref to Candidate object to be recorded in filter object
0093   RecoChargedCandidateRef refMu1;
0094   RecoChargedCandidateRef refMu2;
0095   RecoChargedCandidateRef refTrk;
0096
0097   double e1, e2, e3;
0098   Particle::LorentzVector p, p1, p2, p3;
0099
0100   if (mucands->size() < 2)
0101     return;
0102   if (trkcands->empty())
0103     return;
0104
0105   RecoChargedCandidateCollection::const_iterator mucand1;
0106   RecoChargedCandidateCollection::const_iterator mucand2;
0107   RecoChargedCandidateCollection::const_iterator trkcand;
0108
0109   // get the objects passing the previous filter
0110   Handle<TriggerFilterObjectWithRefs> previousCands;
0111   iEvent.getByToken(previousCandToken_, previousCands);
0112
0113   vector<RecoChargedCandidateRef> vPrevCands;
0114   previousCands->getObjects(TriggerMuon, vPrevCands);
0115
0116   for (mucand1 = mucands->begin(); mucand1 != mucands->end(); ++mucand1) {
0117     TrackRef trk1 = mucand1->get<TrackRef>();
0118     LogDebug("HLTmumutkVtxProducer") << " 1st muon: q*pt= " << trk1->charge() * trk1->pt() << ", eta= " << trk1->eta()
0119                                      << ", hits= " << trk1->numberOfValidHits();
0120
0121     //first check if this muon passed the previous filter
0122     if (!checkPreviousCand(trk1, vPrevCands))
0123       continue;
0124
0125     // eta and pt cut
0126     if (fabs(trk1->eta()) > maxEta_)
0127       continue;
0128     if (trk1->pt() < minPt_)
0129       continue;
0130
0131     mucand2 = mucand1;
0132     ++mucand2;
0133     for (; mucand2 != mucands->end(); mucand2++) {
0134       TrackRef trk2 = mucand2->get<TrackRef>();
0135
0136       LogDebug("HLTDisplacedMumukFilter") << " 2nd muon: q*pt= " << trk2->charge() * trk2->pt()
0137                                           << ", eta= " << trk2->eta() << ", hits= " << trk2->numberOfValidHits();
0138
0139       //first check if this muon passed the previous filter
0140       if (!checkPreviousCand(trk2, vPrevCands))
0141         continue;
0142       // eta and pt cut
0143       if (fabs(trk2->eta()) > maxEta_)
0144         continue;
0145       if (trk2->pt() < minPt_)
0146         continue;
0147
0148       //loop on track collection
0149       for (trkcand = trkcands->begin(); trkcand != trkcands->end(); ++trkcand) {
0150         TrackRef trk3 = trkcand->get<TrackRef>();
0151         if (overlap(trk1, trk3))
0152           continue;
0153         if (overlap(trk2, trk3))
0154           continue;
0155
0156         LogDebug("HLTDisplacedMumukFilter") << " 3rd track: q*pt= " << trk3->charge() * trk3->pt()
0157                                             << ", eta= " << trk3->eta() << ", hits= " << trk3->numberOfValidHits();
0158
0159         // eta and pt cut
0160         if (fabs(trk3->eta()) > maxEta_)
0161           continue;
0162         if (trk3->pt() < minPt_)
0163           continue;
0164
0165         // Combined system
0166         e1 = sqrt(trk1->momentum().Mag2() + MuMass2);
0167         e2 = sqrt(trk2->momentum().Mag2() + MuMass2);
0168         e3 = sqrt(trk3->momentum().Mag2() + thirdTrackMass2);
0169
0170         p1 = Particle::LorentzVector(trk1->px(), trk1->py(), trk1->pz(), e1);
0171         p2 = Particle::LorentzVector(trk2->px(), trk2->py(), trk2->pz(), e2);
0172         p3 = Particle::LorentzVector(trk3->px(), trk3->py(), trk3->pz(), e3);
0173
0174         p = p1 + p2 + p3;
0175
0176         //invariant mass cut
0177         double invmass = abs(p.mass());
0178         LogDebug("HLTDisplacedMumukFilter") << " Invmass= " << invmass;
0179         if (invmass < minInvMass_)
0180           continue;
0181         if (invmass > maxInvMass_)
0182           continue;
0183
0184         // do the vertex fit
0185         vector<TransientTrack> t_tks;
0186         t_tks.push_back((*theB).build(&trk1));
0187         t_tks.push_back((*theB).build(&trk2));
0188         t_tks.push_back((*theB).build(&trk3));
0189         if (t_tks.size() != 3)
0190           continue;
0191
0192         FreeTrajectoryState InitialFTS = initialFreeState(*trk3, magField);
0193         TrajectoryStateClosestToBeamLine tscb(blsBuilder(InitialFTS, *recoBeamSpotHandle));
0194         double d0sig = tscb.transverseImpactParameter().significance();
0195         if (d0sig < minD0Significance_)
0196           continue;
0197
0198         KalmanVertexFitter kvf;
0199         TransientVertex tv = kvf.vertex(t_tks);
0200         if (!tv.isValid())
0201           continue;
0202         Vertex vertex = tv;
0203
0204         // put vertex in the event
0205         vertexCollection->push_back(vertex);
0206       }
0207     }
0208   }
0209   iEvent.put(std::move(vertexCollection));
0210 }
0211
0212 FreeTrajectoryState HLTmumutkVtxProducer::initialFreeState(const reco::Track& tk, const MagneticField* field) {
0213   Basic3DVector<float> pos(tk.vertex());
0214   GlobalPoint gpos(pos);
0215   Basic3DVector<float> mom(tk.momentum());
0216   GlobalVector gmom(mom);
0217   GlobalTrajectoryParameters par(gpos, gmom, tk.charge(), field);
0218   CurvilinearTrajectoryError err(tk.covariance());
0219   return FreeTrajectoryState(par, err);
0220 }
0221
0222 bool HLTmumutkVtxProducer::overlap(const TrackRef& trackref1, const TrackRef& trackref2) {
0223   return (reco::deltaR2(trackref1->eta(), trackref1->phi(), trackref2->eta(), trackref2->phi()) < overlapDR2_);
0224 }
0225
0226 bool HLTmumutkVtxProducer::checkPreviousCand(const TrackRef& trackref,
0227                                              const vector<RecoChargedCandidateRef>& refVect) const {
0228   bool ok = false;
0229   for (auto& i : refVect) {
0230     if (i->get<TrackRef>() == trackref) {
0231       ok = true;
0232       break;
0233     }
0234   }
0235   return ok;
0236 }