btau/plugins/HLTDisplacedtktktkVtxProducer.cc

0001 #include <iostream>
0002
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0006 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0007 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0008 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0009 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0010 #include "DataFormats/Candidate/interface/Candidate.h"
0011 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0012 #include "DataFormats/HLTReco/interface/TriggerRefsCollections.h"
0013 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0017 #include "DataFormats/Common/interface/RefToBase.h"
0018 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
0019 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0020 #include "HLTDisplacedtktktkVtxProducer.h"
0021
0022 using namespace edm;
0023 using namespace reco;
0024 using namespace std;
0025 using namespace trigger;
0026 //
0027 // constructors and destructor
0028 //
0029 HLTDisplacedtktktkVtxProducer::HLTDisplacedtktktkVtxProducer(const edm::ParameterSet& iConfig)
0030     : transientTrackRecordToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0031       srcTag_(iConfig.getParameter<edm::InputTag>("Src")),
0032       srcToken_(consumes<reco::RecoChargedCandidateCollection>(srcTag_)),
0033       previousCandTag_(iConfig.getParameter<edm::InputTag>("PreviousCandTag")),
0034       previousCandToken_(consumes<trigger::TriggerFilterObjectWithRefs>(previousCandTag_)),
0035       maxEta_(iConfig.getParameter<double>("MaxEtaTk")),
0036       minPtTk1_(iConfig.getParameter<double>("MinPtResTk1")),
0037       minPtTk2_(iConfig.getParameter<double>("MinPtResTk2")),
0038       minPtTk3_(iConfig.getParameter<double>("MinPtThirdTk")),
0039       minPtRes_(iConfig.getParameter<double>("MinPtRes")),
0040       minPtTri_(iConfig.getParameter<double>("MinPtTri")),
0041       minInvMassRes_(iConfig.getParameter<double>("MinInvMassRes")),
0042       maxInvMassRes_(iConfig.getParameter<double>("MaxInvMassRes")),
0043       minInvMass_(iConfig.getParameter<double>("MinInvMass")),
0044       maxInvMass_(iConfig.getParameter<double>("MaxInvMass")),
0045       massParticle1_(iConfig.getParameter<double>("massParticleRes1")),
0046       massParticle2_(iConfig.getParameter<double>("massParticleRes2")),
0047       massParticle3_(iConfig.getParameter<double>("massParticle3")),
0048       chargeOpt_(iConfig.getParameter<int>("ChargeOpt")),
0049       resOpt_(iConfig.getParameter<int>("ResOpt")),
0050       triggerTypeDaughters_(iConfig.getParameter<int>("triggerTypeDaughters")) {
0051   produces<VertexCollection>();
0052
0053   firstTrackMass = massParticle1_;
0054   secondTrackMass = massParticle2_;
0055   thirdTrackMass = massParticle3_;
0056   firstTrackPt = minPtTk1_;
0057   secondTrackPt = minPtTk2_;
0058   thirdTrackPt = minPtTk3_;
0059   if (resOpt_ <= 0 && massParticle1_ != massParticle2_) {
0060     if (massParticle1_ == massParticle3_) {
0061       std::swap(secondTrackMass, thirdTrackMass);
0062       std::swap(secondTrackPt, thirdTrackPt);
0063     }
0064     if (massParticle2_ == massParticle3_) {
0065       std::swap(firstTrackMass, thirdTrackMass);
0066       std::swap(firstTrackPt, thirdTrackPt);
0067     }
0068   }
0069   firstTrackMass2 = firstTrackMass * firstTrackMass;
0070   secondTrackMass2 = secondTrackMass * secondTrackMass;
0071   thirdTrackMass2 = thirdTrackMass * thirdTrackMass;
0072 }
0073
0074 HLTDisplacedtktktkVtxProducer::~HLTDisplacedtktktkVtxProducer() = default;
0075
0076 void HLTDisplacedtktktkVtxProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077   edm::ParameterSetDescription desc;
0078   desc.add<edm::InputTag>("Src", edm::InputTag("hltL3MuonCandidates"));
0079   desc.add<edm::InputTag>("PreviousCandTag", edm::InputTag(""));
0080   desc.add<double>("MaxEtaTk", 2.5);
0081   desc.add<double>("MinPtResTk1", 0.0);
0082   desc.add<double>("MinPtResTk2", 0.0);
0083   desc.add<double>("MinPtThirdTk", 0.0);
0084   desc.add<double>("MinPtRes", 0.0);
0085   desc.add<double>("MinPtTri", 0.0);
0086   desc.add<double>("MinInvMassRes", 1.0);
0087   desc.add<double>("MaxInvMassRes", 20.0);
0088   desc.add<double>("MinInvMass", 1.0);
0089   desc.add<double>("MaxInvMass", 20.0);
0090   desc.add<double>("massParticleRes1", 0.4937);
0091   desc.add<double>("massParticleRes2", 0.4937);
0092   desc.add<double>("massParticle3", 0.1396);
0093   desc.add<int>("ChargeOpt", -1);
0094   desc.add<int>("ResOpt", 1);
0095   desc.add<int>("triggerTypeDaughters", 0);
0096
0097   descriptions.add("hltDisplacedtktktkVtxProducer", desc);
0098 }
0099
0100 // ------------ method called on each new Event  ------------
0101 void HLTDisplacedtktktkVtxProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0102   // get hold of track trks
0103   Handle<RecoChargedCandidateCollection> trackcands;
0104   iEvent.getByToken(srcToken_, trackcands);
0105   if (trackcands->size() < 3)
0106     return;
0107
0108   // get the transient track builder
0109   auto const& theB = iSetup.getHandle(transientTrackRecordToken_);
0110
0111   std::unique_ptr<VertexCollection> vertexCollection(new VertexCollection());
0112
0113   // look at all trackcands,  check cuts and make vertices
0114   double e1, e2, e3;
0115   Particle::LorentzVector p, pres;
0116
0117   RecoChargedCandidateCollection::const_iterator cand1;
0118   RecoChargedCandidateCollection::const_iterator cand2;
0119   RecoChargedCandidateCollection::const_iterator cand3;
0120
0121   // get the objects passing the previous filter
0122   Handle<TriggerFilterObjectWithRefs> previousCands;
0123   iEvent.getByToken(previousCandToken_, previousCands);
0124
0125   vector<RecoChargedCandidateRef> vPrevCands;
0126   previousCands->getObjects(triggerTypeDaughters_, vPrevCands);
0127
0128   std::vector<bool> candComp;
0129   for (cand1 = trackcands->begin(); cand1 != trackcands->end(); cand1++)
0130     candComp.push_back(checkPreviousCand(cand1->get<TrackRef>(), vPrevCands));
0131
0132   for (cand1 = trackcands->begin(); cand1 != trackcands->end(); cand1++) {
0133     TrackRef tk1 = cand1->get<TrackRef>();
0134     LogDebug("HLTDisplacedtktktkVtxProducer") << " 1st track in loop: q*pt= " << cand1->charge() * cand1->pt()
0135                                               << ", eta= " << cand1->eta() << ", hits= " << tk1->numberOfValidHits();
0136
0137     //first check if this track passed the previous filter
0138     if (!candComp[cand1 - trackcands->begin()])
0139       continue;
0140     // if( ! checkPreviousCand( tk1, vPrevCands) ) continue;
0141
0142     // cuts
0143     if (std::abs(cand1->eta()) > maxEta_)
0144       continue;
0145     if (cand1->pt() < firstTrackPt)
0146       continue;
0147
0148     cand2 = trackcands->begin();
0149     if (firstTrackMass == secondTrackMass) {
0150       cand2 = cand1 + 1;
0151     }
0152
0153     for (; cand2 != trackcands->end(); cand2++) {
0154       TrackRef tk2 = cand2->get<TrackRef>();
0155       if (tk1 == tk2)
0156         continue;
0157       LogDebug("HLTDisplacedtktktkVtxProducer")
0158           << " 2nd track in loop: q*pt= " << cand2->charge() * cand2->pt() << ", eta= " << cand2->eta()
0159           << ", hits= " << tk2->numberOfValidHits() << ", d0= " << tk2->d0();
0160
0161       //first check if this track passed the previous filter
0162       if (!candComp[cand2 - trackcands->begin()])
0163         continue;
0164       // if( ! checkPreviousCand( tk2, vPrevCands) ) continue;
0165
0166       // cuts
0167       if (std::abs(cand2->eta()) > maxEta_)
0168         continue;
0169       if (cand2->pt() < secondTrackPt)
0170         continue;
0171
0172       // opposite sign or same sign for resonance
0173       if (resOpt_ > 0) {
0174         if (chargeOpt_ < 0) {
0175           if (cand1->charge() * cand2->charge() > 0)
0176             continue;
0177         } else if (chargeOpt_ > 0) {
0178           if (cand1->charge() * cand2->charge() < 0)
0179             continue;
0180         }
0181       }
0182
0183       //
0184       // Combined ditrack system
0185       e1 = sqrt(cand1->momentum().Mag2() + firstTrackMass2);
0186       e2 = sqrt(cand2->momentum().Mag2() + secondTrackMass2);
0187       pres = Particle::LorentzVector(cand1->px(), cand1->py(), cand1->pz(), e1) +
0188              Particle::LorentzVector(cand2->px(), cand2->py(), cand2->pz(), e2);
0189
0190       if (resOpt_ > 0) {
0191         if (pres.pt() < minPtRes_)
0192           continue;
0193         double invmassRes = std::abs(pres.mass());
0194         LogDebug("HLTDisplacedtktktkVtxProducer") << " ... 1-2 invmass= " << invmassRes;
0195         if (invmassRes < minInvMassRes_)
0196           continue;
0197         if (invmassRes > maxInvMassRes_)
0198           continue;
0199       }
0200
0201       cand3 = trackcands->begin();
0202       if (firstTrackMass == secondTrackMass && firstTrackMass == thirdTrackMass && resOpt_ <= 0) {
0203         cand3 = cand2 + 1;
0204       }
0205
0206       for (; cand3 != trackcands->end(); cand3++) {
0207         TrackRef tk3 = cand3->get<TrackRef>();
0208         if (tk1 == tk3 || tk2 == tk3)
0209           continue;
0210         LogDebug("HLTDisplacedtktktkVtxProducer")
0211             << " 3rd track in loop: q*pt= " << cand3->charge() * cand3->pt() << ", eta= " << cand3->eta()
0212             << ", hits= " << tk3->numberOfValidHits();
0213
0214         //first check if this track passed the previous filter
0215         if (!candComp[cand3 - trackcands->begin()])
0216           continue;
0217         // if( ! checkPreviousCand( tk3, vPrevCands) ) continue;
0218
0219         // cuts
0220         if (std::abs(cand3->eta()) > maxEta_)
0221           continue;
0222         if (cand3->pt() < thirdTrackPt)
0223           continue;
0224
0225         e3 = sqrt(cand3->momentum().Mag2() + thirdTrackMass2);
0226         p = Particle::LorentzVector(cand1->px(), cand1->py(), cand1->pz(), e1) +
0227             Particle::LorentzVector(cand2->px(), cand2->py(), cand2->pz(), e2) +
0228             Particle::LorentzVector(cand3->px(), cand3->py(), cand3->pz(), e3);
0229
0230         if (p.pt() < minPtTri_)
0231           continue;
0232         double invmass = std::abs(p.mass());
0233         LogDebug("HLTDisplacedtktktkVtxProducer") << " ... 1-2-3 invmass= " << invmass;
0234         if (invmass < minInvMass_)
0235           continue;
0236         if (invmass > maxInvMass_)
0237           continue;
0238
0239         // do the vertex fit
0240         vector<TransientTrack> t_tks;
0241         TransientTrack ttkp1 = (*theB).build(&tk1);
0242         TransientTrack ttkp2 = (*theB).build(&tk2);
0243         TransientTrack ttkp3 = (*theB).build(&tk3);
0244
0245         t_tks.push_back(ttkp1);
0246         t_tks.push_back(ttkp2);
0247         t_tks.push_back(ttkp3);
0248
0249         if (t_tks.size() != 3)
0250           continue;
0251
0252         KalmanVertexFitter kvf;
0253         TransientVertex tv = kvf.vertex(t_tks);
0254
0255         if (!tv.isValid())
0256           continue;
0257
0258         Vertex vertex = tv;
0259
0260         // put vertex in the event
0261         vertexCollection->push_back(vertex);
0262       }
0263     }
0264   }
0265   iEvent.put(std::move(vertexCollection));
0266 }
0267
0268 bool HLTDisplacedtktktkVtxProducer::checkPreviousCand(const TrackRef& trackref,
0269                                                       const vector<RecoChargedCandidateRef>& refVect) const {
0270   bool ok = false;
0271   for (auto& i : refVect) {
0272     if (i->get<TrackRef>() == trackref) {
0273       ok = true;
0274       break;
0275     }
0276   }
0277   return ok;
0278 }