Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:02 

0001 #ifndef InclusiveVertexFinder_h
0002 #define InclusiveVertexFinder_h
0003 #include <memory>
0004 
0005 #include "FWCore/Framework/interface/stream/EDProducer.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010 #include "FWCore/Utilities/interface/ESGetToken.h"
0011 
0012 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0013 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0014 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0015 
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019 #include "DataFormats/VertexReco/interface/Vertex.h"
0020 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0021 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0022 
0023 #include "RecoVertex/ConfigurableVertexReco/interface/ConfigurableVertexReconstructor.h"
0024 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
0025 #include "TrackingTools/IPTools/interface/IPTools.h"
0026 #include "RecoVertex/AdaptiveVertexFinder/interface/TracksClusteringFromDisplacedSeed.h"
0027 
0028 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
0029 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h"
0030 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackCompatibilityEstimator.h"
0031 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
0032 #include "RecoVertex/MultiVertexFit/interface/MultiVertexFitter.h"
0033 
0034 #include "RecoVertex/AdaptiveVertexFinder/interface/TTHelpers.h"
0035 #include "RecoVertex/AdaptiveVertexFinder/interface/SVTimeHelpers.h"
0036 #include "FWCore/Utilities/interface/isFinite.h"
0037 
0038 #include <type_traits>
0039 
0040 //#define VTXDEBUG 1
0041 template <class InputContainer, class VTX>
0042 class TemplatedInclusiveVertexFinder : public edm::stream::EDProducer<> {
0043 public:
0044   typedef std::vector<VTX> Product;
0045   typedef typename InputContainer::value_type TRK;
0046   TemplatedInclusiveVertexFinder(const edm::ParameterSet &params);
0047 
0048   static void fillDescriptions(edm::ConfigurationDescriptions &cdesc) {
0049     edm::ParameterSetDescription pdesc;
0050     pdesc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));
0051     pdesc.add<edm::InputTag>("primaryVertices", edm::InputTag("offlinePrimaryVertices"));
0052     if (std::is_same<VTX, reco::Vertex>::value) {
0053       pdesc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0054       pdesc.add<unsigned int>("minHits", 8);
0055     } else if (std::is_same<VTX, reco::VertexCompositePtrCandidate>::value) {
0056       pdesc.add<edm::InputTag>("tracks", edm::InputTag("particleFlow"));
0057       pdesc.add<unsigned int>("minHits", 0);
0058     } else {
0059       pdesc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0060     }
0061 
0062     pdesc.add<double>("maximumLongitudinalImpactParameter", 0.3);
0063     pdesc.add<double>("maximumTimeSignificance", 3.0);
0064     pdesc.add<double>("minPt", 0.8);
0065     pdesc.add<unsigned int>("maxNTracks", 30);
0066     //clusterizer pset
0067     edm::ParameterSetDescription clusterizer;
0068     clusterizer.add<double>("seedMax3DIPSignificance", 9999.0);
0069     clusterizer.add<double>("seedMax3DIPValue", 9999.0);
0070     clusterizer.add<double>("seedMin3DIPSignificance", 1.2);
0071     clusterizer.add<double>("seedMin3DIPValue", 0.005);
0072     clusterizer.add<double>("clusterMaxDistance", 0.05);
0073     clusterizer.add<double>("clusterMaxSignificance", 4.5);
0074     clusterizer.add<double>("distanceRatio", 20.0);
0075     clusterizer.add<double>("clusterMinAngleCosine", 0.5);
0076     clusterizer.add<double>("maxTimeSignificance", 3.5);
0077     pdesc.add<edm::ParameterSetDescription>("clusterizer", clusterizer);
0078     // vertex and fitter config
0079     pdesc.add<double>("vertexMinAngleCosine", 0.95);
0080     pdesc.add<double>("vertexMinDLen2DSig", 2.5);
0081     pdesc.add<double>("vertexMinDLenSig", 0.5);
0082     pdesc.add<double>("fitterSigmacut", 3.0);
0083     pdesc.add<double>("fitterTini", 256.0);
0084     pdesc.add<double>("fitterRatio", 0.25);
0085     pdesc.add<bool>("useDirectVertexFitter", true);
0086     pdesc.add<bool>("useVertexReco", true);
0087     // vertexReco pset
0088     edm::ParameterSetDescription vertexReco;
0089     vertexReco.add<std::string>("finder", std::string("avr"));
0090     vertexReco.add<double>("primcut", 1.0);
0091     vertexReco.add<double>("seccut", 3.0);
0092     vertexReco.add<bool>("smoothing", true);
0093     pdesc.add<edm::ParameterSetDescription>("vertexReco", vertexReco);
0094     if (std::is_same<VTX, reco::Vertex>::value) {
0095       cdesc.add("inclusiveVertexFinderDefault", pdesc);
0096     } else if (std::is_same<VTX, reco::VertexCompositePtrCandidate>::value) {
0097       cdesc.add("inclusiveCandidateVertexFinderDefault", pdesc);
0098     } else {
0099       cdesc.addDefault(pdesc);
0100     }
0101   }
0102 
0103   void produce(edm::Event &event, const edm::EventSetup &es) override;
0104 
0105 private:
0106   bool trackFilter(const reco::Track &track) const;
0107   std::pair<std::vector<reco::TransientTrack>, GlobalPoint> nearTracks(const reco::TransientTrack &seed,
0108                                                                        const std::vector<reco::TransientTrack> &tracks,
0109                                                                        const reco::Vertex &primaryVertex) const;
0110 
0111   edm::EDGetTokenT<reco::BeamSpot> token_beamSpot;
0112   edm::EDGetTokenT<reco::VertexCollection> token_primaryVertex;
0113   edm::EDGetTokenT<InputContainer> token_tracks;
0114   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> token_trackBuilder;
0115   unsigned int minHits;
0116   unsigned int maxNTracks;
0117   double maxLIP;
0118   double maxTimeSig;
0119   double minPt;
0120   double vertexMinAngleCosine;
0121   double vertexMinDLen2DSig;
0122   double vertexMinDLenSig;
0123   double fitterSigmacut;
0124   double fitterTini;
0125   double fitterRatio;
0126   bool useVertexFitter;
0127   bool useVertexReco;
0128   std::unique_ptr<VertexReconstructor> vtxReco;
0129   std::unique_ptr<TracksClusteringFromDisplacedSeed> clusterizer;
0130 };
0131 template <class InputContainer, class VTX>
0132 TemplatedInclusiveVertexFinder<InputContainer, VTX>::TemplatedInclusiveVertexFinder(const edm::ParameterSet &params)
0133     : minHits(params.getParameter<unsigned int>("minHits")),
0134       maxNTracks(params.getParameter<unsigned int>("maxNTracks")),
0135       maxLIP(params.getParameter<double>("maximumLongitudinalImpactParameter")),
0136       maxTimeSig(params.getParameter<double>("maximumTimeSignificance")),
0137       minPt(params.getParameter<double>("minPt")),                                //0.8
0138       vertexMinAngleCosine(params.getParameter<double>("vertexMinAngleCosine")),  //0.98
0139       vertexMinDLen2DSig(params.getParameter<double>("vertexMinDLen2DSig")),      //2.5
0140       vertexMinDLenSig(params.getParameter<double>("vertexMinDLenSig")),          //0.5
0141       fitterSigmacut(params.getParameter<double>("fitterSigmacut")),
0142       fitterTini(params.getParameter<double>("fitterTini")),
0143       fitterRatio(params.getParameter<double>("fitterRatio")),
0144       useVertexFitter(params.getParameter<bool>("useDirectVertexFitter")),
0145       useVertexReco(params.getParameter<bool>("useVertexReco")),
0146       vtxReco(new ConfigurableVertexReconstructor(params.getParameter<edm::ParameterSet>("vertexReco"))),
0147       clusterizer(new TracksClusteringFromDisplacedSeed(params.getParameter<edm::ParameterSet>("clusterizer")))
0148 
0149 {
0150   token_beamSpot = consumes<reco::BeamSpot>(params.getParameter<edm::InputTag>("beamSpot"));
0151   token_primaryVertex = consumes<reco::VertexCollection>(params.getParameter<edm::InputTag>("primaryVertices"));
0152   token_tracks = consumes<InputContainer>(params.getParameter<edm::InputTag>("tracks"));
0153   token_trackBuilder =
0154       esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0155   produces<Product>();
0156   //produces<reco::VertexCollection>("multi");
0157 }
0158 template <class InputContainer, class VTX>
0159 bool TemplatedInclusiveVertexFinder<InputContainer, VTX>::trackFilter(const reco::Track &track) const {
0160   if (track.hitPattern().numberOfValidHits() < (int)minHits)
0161     return false;
0162   if (track.pt() < minPt)
0163     return false;
0164 
0165   return true;
0166 }
0167 
0168 template <class InputContainer, class VTX>
0169 void TemplatedInclusiveVertexFinder<InputContainer, VTX>::produce(edm::Event &event, const edm::EventSetup &es) {
0170   using namespace reco;
0171 
0172   VertexDistance3D vdist;
0173   VertexDistanceXY vdist2d;
0174   MultiVertexFitter theMultiVertexFitter;
0175   AdaptiveVertexFitter theAdaptiveFitter(GeometricAnnealing(fitterSigmacut, fitterTini, fitterRatio),
0176                                          DefaultLinearizationPointFinder(),
0177                                          KalmanVertexUpdator<5>(),
0178                                          KalmanVertexTrackCompatibilityEstimator<5>(),
0179                                          KalmanVertexSmoother());
0180 
0181   edm::Handle<BeamSpot> beamSpot;
0182   event.getByToken(token_beamSpot, beamSpot);
0183 
0184   edm::Handle<VertexCollection> primaryVertices;
0185   event.getByToken(token_primaryVertex, primaryVertices);
0186 
0187   edm::Handle<InputContainer> tracks;
0188   event.getByToken(token_tracks, tracks);
0189 
0190   edm::ESHandle<TransientTrackBuilder> trackBuilder = es.getHandle(token_trackBuilder);
0191 
0192   auto recoVertices = std::make_unique<Product>();
0193   if (!primaryVertices->empty()) {
0194     const reco::Vertex &pv = (*primaryVertices)[0];
0195     GlobalPoint ppv(pv.position().x(), pv.position().y(), pv.position().z());
0196 
0197     std::vector<TransientTrack> tts;
0198     //Fill transient track vector
0199     for (typename InputContainer::const_iterator track = tracks->begin(); track != tracks->end(); ++track) {
0200       //TransientTrack tt = trackBuilder->build(ref);
0201       //TrackRef ref(tracks, track - tracks->begin());
0202       TransientTrack tt(tthelpers::buildTT(tracks, trackBuilder, track - tracks->begin()));
0203       if (!tt.isValid())
0204         continue;
0205       if (!trackFilter(tt.track()))
0206         continue;
0207       if (std::abs(tt.track().dz(pv.position())) > maxLIP)
0208         continue;
0209       if (edm::isFinite(tt.timeExt()) && pv.covariance(3, 3) > 0.) {  // only apply if time available
0210         auto tError = std::sqrt(std::pow(tt.dtErrorExt(), 2) + pv.covariance(3, 3));
0211         auto dtSig = std::abs(tt.timeExt() - pv.t()) / tError;
0212         if (dtSig > maxTimeSig)
0213           continue;
0214       }
0215       tt.setBeamSpot(*beamSpot);
0216       tts.push_back(tt);
0217     }
0218     std::vector<TracksClusteringFromDisplacedSeed::Cluster> clusters = clusterizer->clusters(pv, tts);
0219 
0220     //Create BS object from PV to feed in the AVR
0221     BeamSpot::CovarianceMatrix cov;
0222     for (unsigned int i = 0; i < 7; i++) {
0223       for (unsigned int j = 0; j < 7; j++) {
0224         if (i < 3 && j < 3)
0225           cov(i, j) = pv.covariance(i, j);
0226         else
0227           cov(i, j) = 0.0;
0228       }
0229     }
0230     BeamSpot bs(pv.position(), 0.0, 0.0, 0.0, 0.0, cov, BeamSpot::Unknown);
0231 
0232     int i = 0;
0233 #ifdef VTXDEBUG
0234 
0235     std::cout << "CLUSTERS " << clusters.size() << std::endl;
0236 #endif
0237 
0238     for (std::vector<TracksClusteringFromDisplacedSeed::Cluster>::iterator cluster = clusters.begin();
0239          cluster != clusters.end();
0240          ++cluster, ++i) {
0241       if (cluster->tracks.size() < 2 || cluster->tracks.size() > maxNTracks)
0242         continue;
0243       std::vector<TransientVertex> vertices;
0244       if (useVertexReco) {
0245         vertices = vtxReco->vertices(cluster->tracks, bs);  // attempt with config given reconstructor
0246       }
0247       TransientVertex singleFitVertex;
0248       if (useVertexFitter) {
0249         singleFitVertex = theAdaptiveFitter.vertex(cluster->tracks, cluster->seedPoint);  //attempt with direct fitting
0250         if (singleFitVertex.isValid())
0251           vertices.push_back(singleFitVertex);
0252       }
0253 
0254       // for each transient vertex state determine if a time can be measured and fill covariance
0255       if (pv.covariance(3, 3) > 0.) {
0256         for (auto &vtx : vertices) {
0257           svhelper::updateVertexTime(vtx);
0258         }
0259       }
0260 
0261       for (std::vector<TransientVertex>::const_iterator v = vertices.begin(); v != vertices.end(); ++v) {
0262         Measurement1D dlen = vdist.distance(pv, *v);
0263         Measurement1D dlen2 = vdist2d.distance(pv, *v);
0264 #ifdef VTXDEBUG
0265         VTX vv(*v);
0266         std::cout << "V chi2/n: " << v->normalisedChiSquared() << " ndof: " << v->degreesOfFreedom();
0267         std::cout << " dlen: " << dlen.value() << " error: " << dlen.error() << " signif: " << dlen.significance();
0268         std::cout << " dlen2: " << dlen2.value() << " error2: " << dlen2.error()
0269                   << " signif2: " << dlen2.significance();
0270         std::cout << " pos: " << vv.position() << " error: " << vv.xError() << " " << vv.yError() << " " << vv.zError()
0271                   << std::endl;
0272         std::cout << " time: " << vv.time() << " error: " << vv.tError() << std::endl;
0273 #endif
0274         GlobalVector dir;
0275         std::vector<reco::TransientTrack> ts = v->originalTracks();
0276         for (std::vector<reco::TransientTrack>::const_iterator i = ts.begin(); i != ts.end(); ++i) {
0277           float w = v->trackWeight(*i);
0278           if (w > 0.5)
0279             dir += i->impactPointState().globalDirection();
0280 #ifdef VTXDEBUG
0281           std::cout << "\t[" << (*i).track().pt() << ": " << (*i).track().eta() << ", " << (*i).track().phi() << "], "
0282                     << w << std::endl;
0283 #endif
0284         }
0285         GlobalPoint sv((*v).position().x(), (*v).position().y(), (*v).position().z());
0286         float vscal = dir.unit().dot((sv - ppv).unit());
0287         if (dlen.significance() > vertexMinDLenSig &&
0288             ((vertexMinAngleCosine > 0) ? (vscal > vertexMinAngleCosine) : (vscal < vertexMinAngleCosine)) &&
0289             v->normalisedChiSquared() < 10 && dlen2.significance() > vertexMinDLen2DSig) {
0290           recoVertices->push_back(*v);
0291 
0292 #ifdef VTXDEBUG
0293           std::cout << "ADDED" << std::endl;
0294 #endif
0295         }
0296       }
0297     }
0298 #ifdef VTXDEBUG
0299 
0300     std::cout << "Final put  " << recoVertices->size() << std::endl;
0301 #endif
0302   }
0303 
0304   event.put(std::move(recoVertices));
0305 }
0306 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!