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