Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:38:05 

0001 // -*- C++ -*-
0002 //
0003 // Package:    PixelVertexProducer
0004 // Class:      PixelVertexProducer
0005 //
0006 /**\class PixelVertexProducer PixelVertexProducer.h PixelVertexFinding/interface/PixelVertexProducer.h
0007 
0008  Description: This produces 1D (z only) primary vertexes using only pixel information.
0009 
0010  Implementation:
0011      This producer can use either the Divisive Primary Vertex Finder
0012      or the Histogramming Primary Vertex Finder (currently not
0013      implemented).  It relies on the PixelTripletProducer and
0014      PixelTrackProducer having already been run upstream.   This is
0015      code ported from ORCA originally written by S Cucciarelli, M
0016      Konecki, D Kotlinski.
0017 */
0018 //
0019 // Original Author:  Aaron Dominguez (UNL)
0020 //         Created:  Thu May 25 10:17:32 CDT 2006
0021 //
0022 
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/Utilities/interface/InputTag.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0032 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0033 #include "DataFormats/VertexReco/interface/Vertex.h"
0034 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0035 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0036 #include "RecoVertex/PixelVertexFinding/interface/DivisiveVertexFinder.h"
0037 #include <memory>
0038 #include <string>
0039 #include <cmath>
0040 
0041 class PixelVertexProducer : public edm::stream::EDProducer<> {
0042 public:
0043   explicit PixelVertexProducer(const edm::ParameterSet&);
0044   ~PixelVertexProducer() override;
0045 
0046   void produce(edm::Event&, const edm::EventSetup&) override;
0047 
0048 private:
0049   // ----------member data ---------------------------
0050   // Turn on debug printing if verbose_ > 0
0051   const int verbose_;
0052   // Tracking cuts before sending tracks to vertex algo
0053   const double ptMin_;
0054   const bool method2;
0055   const edm::InputTag trackCollName;
0056   const edm::EDGetTokenT<reco::TrackCollection> token_Tracks;
0057   const edm::EDGetTokenT<reco::BeamSpot> token_BeamSpot;
0058 
0059   DivisiveVertexFinder* dvf_;
0060 };
0061 
0062 PixelVertexProducer::PixelVertexProducer(const edm::ParameterSet& conf)
0063     // 0 silent, 1 chatty, 2 loud
0064     : verbose_(conf.getParameter<int>("Verbosity")),
0065       // 1.0 GeV
0066       ptMin_(conf.getParameter<double>("PtMin")),
0067       method2(conf.getParameter<bool>("Method2")),
0068       trackCollName(conf.getParameter<edm::InputTag>("TrackCollection")),
0069       token_Tracks(consumes<reco::TrackCollection>(trackCollName)),
0070       token_BeamSpot(consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpot"))) {
0071   // Register my product
0072   produces<reco::VertexCollection>();
0073 
0074   // Setup shop
0075   std::string finder = conf.getParameter<std::string>("Finder");  // DivisiveVertexFinder
0076   bool useError = conf.getParameter<bool>("UseError");            // true
0077   bool wtAverage = conf.getParameter<bool>("WtAverage");          // true
0078   double zOffset = conf.getParameter<double>("ZOffset");          // 5.0 sigma
0079   double zSeparation = conf.getParameter<double>("ZSeparation");  // 0.05 cm
0080   int ntrkMin = conf.getParameter<int>("NTrkMin");                // 3
0081   // Tracking requirements before sending a track to be considered for vtx
0082 
0083   double track_pt_min = ptMin_;
0084   double track_pt_max = 10.;
0085   double track_chi2_max = 9999999.;
0086   double track_prob_min = -1.;
0087 
0088   if (conf.exists("PVcomparer")) {
0089     edm::ParameterSet PVcomparerPSet = conf.getParameter<edm::ParameterSet>("PVcomparer");
0090     track_pt_min = PVcomparerPSet.getParameter<double>("track_pt_min");
0091     if (track_pt_min != ptMin_) {
0092       if (track_pt_min < ptMin_)
0093         edm::LogInfo("PixelVertexProducer")
0094             << "minimum track pT setting differs between PixelVertexProducer (" << ptMin_ << ") and PVcomparer ("
0095             << track_pt_min << ") [PVcomparer considers tracks w/ lower threshold than PixelVertexProducer does] !!!";
0096       else
0097         edm::LogInfo("PixelVertexProducer") << "minimum track pT setting differs between PixelVertexProducer ("
0098                                             << ptMin_ << ") and PVcomparer (" << track_pt_min << ") !!!";
0099     }
0100     track_pt_max = PVcomparerPSet.getParameter<double>("track_pt_max");
0101     track_chi2_max = PVcomparerPSet.getParameter<double>("track_chi2_max");
0102     track_prob_min = PVcomparerPSet.getParameter<double>("track_prob_min");
0103   }
0104 
0105   if (finder == "DivisiveVertexFinder") {
0106     if (verbose_ > 0)
0107       edm::LogInfo("PixelVertexProducer") << ": Using the DivisiveVertexFinder\n";
0108     dvf_ = new DivisiveVertexFinder(track_pt_min,
0109                                     track_pt_max,
0110                                     track_chi2_max,
0111                                     track_prob_min,
0112                                     zOffset,
0113                                     ntrkMin,
0114                                     useError,
0115                                     zSeparation,
0116                                     wtAverage,
0117                                     verbose_);
0118   } else {  // Finder not supported, or you made a mistake in your request
0119     // throw an exception once I figure out how CMSSW does this
0120   }
0121 }
0122 
0123 PixelVertexProducer::~PixelVertexProducer() { delete dvf_; }
0124 
0125 void PixelVertexProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0126   // First fish the pixel tracks out of the event
0127   edm::Handle<reco::TrackCollection> trackCollection;
0128   e.getByToken(token_Tracks, trackCollection);
0129   const reco::TrackCollection tracks = *(trackCollection.product());
0130   if (verbose_ > 0)
0131     edm::LogInfo("PixelVertexProducer") << ": Found " << tracks.size() << " tracks in TrackCollection called "
0132                                         << trackCollName << "\n";
0133 
0134   // Second, make a collection of pointers to the tracks we want for the vertex finder
0135   reco::TrackRefVector trks;
0136   for (unsigned int i = 0; i < tracks.size(); i++) {
0137     if (tracks[i].pt() > ptMin_)
0138       trks.push_back(reco::TrackRef(trackCollection, i));
0139   }
0140   if (verbose_ > 0)
0141     edm::LogInfo("PixelVertexProducer") << ": Selected " << trks.size() << " of these tracks for vertexing\n";
0142 
0143   edm::Handle<reco::BeamSpot> bsHandle;
0144   e.getByToken(token_BeamSpot, bsHandle);
0145   math::XYZPoint myPoint(0., 0., 0.);
0146   if (bsHandle.isValid())
0147     //FIXME: fix last coordinate with vertex.z() at same time
0148     myPoint = math::XYZPoint(bsHandle->x0(), bsHandle->y0(), 0.);
0149 
0150   // Third, ship these tracks off to be vertexed
0151   auto vertexes = std::make_unique<reco::VertexCollection>();
0152   bool ok;
0153   if (method2) {
0154     ok = dvf_->findVertexesAlt(trks,  // input
0155                                *vertexes,
0156                                myPoint);  // output
0157     if (verbose_ > 0)
0158       edm::LogInfo("PixelVertexProducer") << "Method2 returned status of " << ok;
0159   } else {
0160     ok = dvf_->findVertexes(trks,        // input
0161                             *vertexes);  // output
0162     if (verbose_ > 0)
0163       edm::LogInfo("PixelVertexProducer") << "Method1 returned status of " << ok;
0164   }
0165 
0166   if (verbose_ > 0) {
0167     edm::LogInfo("PixelVertexProducer") << ": Found " << vertexes->size() << " vertexes\n";
0168     for (unsigned int i = 0; i < vertexes->size(); ++i) {
0169       edm::LogInfo("PixelVertexProducer")
0170           << "Vertex number " << i << " has " << (*vertexes)[i].tracksSize() << " tracks with a position of "
0171           << (*vertexes)[i].z() << " +- " << std::sqrt((*vertexes)[i].covariance(2, 2));
0172     }
0173   }
0174 
0175   if (bsHandle.isValid()) {
0176     const reco::BeamSpot& bs = *bsHandle;
0177 
0178     for (unsigned int i = 0; i < vertexes->size(); ++i) {
0179       double z = (*vertexes)[i].z();
0180       double x = bs.x0() + bs.dxdz() * (z - bs.z0());
0181       double y = bs.y0() + bs.dydz() * (z - bs.z0());
0182       reco::Vertex v(reco::Vertex::Point(x, y, z),
0183                      (*vertexes)[i].error(),
0184                      (*vertexes)[i].chi2(),
0185                      (*vertexes)[i].ndof(),
0186                      (*vertexes)[i].tracksSize());
0187       //Copy also the tracks
0188       for (std::vector<reco::TrackBaseRef>::const_iterator it = (*vertexes)[i].tracks_begin();
0189            it != (*vertexes)[i].tracks_end();
0190            it++) {
0191         v.add(*it);
0192       }
0193       (*vertexes)[i] = v;
0194     }
0195   } else {
0196     edm::LogWarning("PixelVertexProducer") << "No beamspot found. Using returning vertexes with (0,0,Z) ";
0197   }
0198 
0199   if (vertexes->empty() && bsHandle.isValid()) {
0200     const reco::BeamSpot& bs = *bsHandle;
0201 
0202     GlobalError bse(bs.rotatedCovariance3D());
0203     if ((bse.cxx() <= 0.) || (bse.cyy() <= 0.) || (bse.czz() <= 0.)) {
0204       AlgebraicSymMatrix33 we;
0205       we(0, 0) = 10000;
0206       we(1, 1) = 10000;
0207       we(2, 2) = 10000;
0208       vertexes->push_back(reco::Vertex(bs.position(), we, 0., 0., 0));
0209 
0210       edm::LogInfo("PixelVertexProducer")
0211           << "No vertices found. Beamspot with invalid errors " << bse.matrix() << std::endl
0212           << "Will put Vertex derived from dummy-fake BeamSpot into Event.\n"
0213           << (*vertexes)[0].x() << "\n"
0214           << (*vertexes)[0].y() << "\n"
0215           << (*vertexes)[0].z() << "\n";
0216     } else {
0217       vertexes->push_back(reco::Vertex(bs.position(), bs.rotatedCovariance3D(), 0., 0., 0));
0218 
0219       edm::LogInfo("PixelVertexProducer") << "No vertices found. Will put Vertex derived from BeamSpot into Event:\n"
0220                                           << (*vertexes)[0].x() << "\n"
0221                                           << (*vertexes)[0].y() << "\n"
0222                                           << (*vertexes)[0].z() << "\n";
0223     }
0224   }
0225 
0226   else if (vertexes->empty() && !bsHandle.isValid()) {
0227     edm::LogWarning("PixelVertexProducer") << "No beamspot and no vertex found. No vertex returned.";
0228   }
0229 
0230   e.put(std::move(vertexes));
0231 }
0232 
0233 DEFINE_FWK_MODULE(PixelVertexProducer);
[ 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!