HiTracking/plugins/HIBestVertexProducer.cc

0001 #include "RecoHI/HiTracking/interface/HIBestVertexProducer.h"
0002
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <algorithm>
0012 #include <iostream>
0013 using namespace std;
0014 using namespace edm;
0015
0016 /*****************************************************************************/
0017 HIBestVertexProducer::HIBestVertexProducer(const edm::ParameterSet& ps)
0018     : theConfig(ps),
0019       theBeamSpotTag(consumes<reco::BeamSpot>(ps.getParameter<edm::InputTag>("beamSpotLabel"))),
0020       theMedianVertexCollection(
0021           consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("medianVertexCollection"))),
0022       theAdaptiveVertexCollection(
0023           consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("adaptiveVertexCollection"))),
0024       theUseFinalAdaptiveVertexCollection(ps.getParameter<bool>("useFinalAdaptiveVertexCollection")) {
0025   if (theUseFinalAdaptiveVertexCollection) {
0026     theFinalAdaptiveVertexCollection =
0027         consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("finalAdaptiveVertexCollection"));
0028   }
0029   produces<reco::VertexCollection>();
0030 }
0031
0032 /*****************************************************************************/
0033 HIBestVertexProducer::~HIBestVertexProducer() {}
0034
0035 /*****************************************************************************/
0036 void HIBestVertexProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0037   edm::ParameterSetDescription desc;
0038   desc.add<InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
0039   desc.add<InputTag>("adaptiveVertexCollection", edm::InputTag("hiBestAdaptiveVertex"));
0040   desc.add<InputTag>("medianVertexCollection", edm::InputTag("hiPixelMedianVertex"));
0041   desc.add<bool>("useFinalAdaptiveVertexCollection", false);
0042   desc.add<InputTag>("finalAdaptiveVertexCollection", edm::InputTag("hiBestOfflinePrimaryVertex"));
0043
0044   descriptions.add("hiSelectedPixelVertex", desc);
0045 }
0046
0047 /*****************************************************************************/
0048 void HIBestVertexProducer::beginJob() {}
0049
0050 /*****************************************************************************/
0051 void HIBestVertexProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
0052   // 1. use best adaptive vertex preferentially
0053   // 2. use median vertex in case where adaptive algorithm failed
0054   // 3. use beamspot if netither vertexing method succeeds
0055
0056   // New vertex collection
0057   auto newVertexCollection = std::make_unique<reco::VertexCollection>();
0058
0059   //** Get precise adaptive vertex **/
0060   edm::Handle<reco::VertexCollection> vc1;
0061   ev.getByToken(theAdaptiveVertexCollection, vc1);
0062   const reco::VertexCollection* vertices1 = vc1.product();
0063
0064   if (vertices1->empty())
0065     LogError("HeavyIonVertexing") << "adaptive vertex collection is empty!" << endl;
0066
0067   //** Final vertex collection if needed **//
0068   const double maxZError = 3.0;  //any vtx with error > this number is considered no good
0069   bool hasFinalVertex = false;
0070   if (theUseFinalAdaptiveVertexCollection) {
0071     edm::Handle<reco::VertexCollection> vc0;
0072     ev.getByToken(theFinalAdaptiveVertexCollection, vc0);
0073     const reco::VertexCollection* vertices0 = vc0.product();
0074     if (vertices0->empty())
0075       LogInfo("HeavyIonVertexing") << "final adaptive vertex collection is empty!" << endl;
0076
0077     //if using final vertex and has a good vertex, use it
0078     if (vertices0->begin()->zError() < maxZError) {
0079       hasFinalVertex = true;
0080       auto const& vertex0 = vertices0->front();
0081       newVertexCollection->push_back(vertex0);
0082       LogInfo("HeavyIonVertexing") << "adaptive vertex:\n vz = (" << vertex0.x() << ", " << vertex0.y() << ", "
0083                                    << vertex0.z() << ")"
0084                                    << "\n error = (" << vertex0.xError() << ", " << vertex0.yError() << ", "
0085                                    << vertex0.zError() << ")" << endl;
0086     }
0087   }
0088
0089   //otherwise use the pixel track adaptive vertex if it is good
0090   if (!hasFinalVertex) {
0091     if (vertices1->begin()->zError() < maxZError) {
0092       reco::VertexCollection::const_iterator vertex1 = vertices1->begin();
0093       newVertexCollection->push_back(*vertex1);
0094
0095       LogInfo("HeavyIonVertexing") << "adaptive vertex:\n vz = (" << vertex1->x() << ", " << vertex1->y() << ", "
0096                                    << vertex1->z() << ")"
0097                                    << "\n error = (" << vertex1->xError() << ", " << vertex1->yError() << ", "
0098                                    << vertex1->zError() << ")" << endl;
0099     } else {
0100       //** Get fast median vertex **/
0101       edm::Handle<reco::VertexCollection> vc2;
0102       ev.getByToken(theMedianVertexCollection, vc2);
0103       const reco::VertexCollection* vertices2 = vc2.product();
0104
0105       //** Get beam spot position and error **/
0106       reco::BeamSpot beamSpot;
0107       edm::Handle<reco::BeamSpot> beamSpotHandle;
0108       ev.getByToken(theBeamSpotTag, beamSpotHandle);
0109
0110       if (beamSpotHandle.isValid())
0111         beamSpot = *beamSpotHandle;
0112       else
0113         LogError("HeavyIonVertexing") << "no beamspot found " << endl;
0114
0115       if (!vertices2->empty()) {
0116         reco::VertexCollection::const_iterator vertex2 = vertices2->begin();
0117         reco::Vertex::Error err;
0118         err(0, 0) = pow(beamSpot.BeamWidthX(), 2);
0119         err(1, 1) = pow(beamSpot.BeamWidthY(), 2);
0120         err(2, 2) = pow(vertex2->zError(), 2);
0121         reco::Vertex newVertex(reco::Vertex::Point(beamSpot.x0(), beamSpot.y0(), vertex2->z()), err, 0, 1, 1);
0122         newVertexCollection->push_back(newVertex);
0123
0124         LogInfo("HeavyIonVertexing") << "median vertex + beamspot: \n position = (" << newVertex.x() << ", "
0125                                      << newVertex.y() << ", " << newVertex.z() << ")"
0126                                      << "\n error = (" << newVertex.xError() << ", " << newVertex.yError() << ", "
0127                                      << newVertex.zError() << ")" << endl;
0128
0129       } else {
0130         reco::Vertex::Error err;
0131         err(0, 0) = pow(beamSpot.BeamWidthX(), 2);
0132         err(1, 1) = pow(beamSpot.BeamWidthY(), 2);
0133         err(2, 2) = pow(beamSpot.sigmaZ(), 2);
0134         reco::Vertex newVertex(beamSpot.position(), err, 0, 0, 1);
0135         newVertexCollection->push_back(newVertex);
0136
0137         LogInfo("HeavyIonVertexing") << "beam spot: \n position = (" << newVertex.x() << ", " << newVertex.y() << ", "
0138                                      << newVertex.z() << ")"
0139                                      << "\n error = (" << newVertex.xError() << ", " << newVertex.yError() << ", "
0140                                      << newVertex.zError() << ")" << endl;
0141       }
0142     }
0143   }
0144
0145   // put new vertex collection into event
0146   ev.put(std::move(newVertexCollection));
0147 }