DivisiveVertexFinder.cc

Go to the documentation of this file.

#include "RecoPixelVertexing/PixelVertexFinding/interface/DivisiveVertexFinder.h"
#include "RecoPixelVertexing/PixelVertexFinding/interface/PVPositionBuilder.h"
#include "RecoPixelVertexing/PixelVertexFinding/interface/PVCluster.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include <utility>
#include <vector>
#include <map>
#include <algorithm>
#include <cmath>

DivisiveVertexFinder::DivisiveVertexFinder(double track_pt_min,
                                           double track_pt_max,
                                           double track_chi2_max,
                                           double track_prob_min,
                                           double zOffset,
                                           int ntrkMin,
                                           bool useError,
                                           double zSeparation,
                                           bool wtAverage,
                                           int verbosity)
    : zOffset_(zOffset),
      zSeparation_(zSeparation),
      ntrkMin_(ntrkMin),
      useError_(useError),
      wtAverage_(wtAverage),
      divmeth_(zOffset, ntrkMin, useError, zSeparation, wtAverage),
      verbose_(verbosity) {
  pvComparer_ = new PVClusterComparer(track_pt_min, track_pt_max, track_chi2_max, track_prob_min);
}

DivisiveVertexFinder::~DivisiveVertexFinder() {}

bool DivisiveVertexFinder::findVertexes(const reco::TrackRefVector &trks,    // input
                                        reco::VertexCollection &vertexes) {  // output
  PVPositionBuilder pos;
  Measurement1D vz;
  if (wtAverage_) {
    vz = pos.wtAverage(trks);
  } else {
    vz = pos.average(trks);
  }
  reco::Vertex::Error err;
  err(2, 2) = vz.error() * vz.error();

  reco::Vertex v(reco::Vertex::Point(0, 0, vz.value()), err, 0, 1, trks.size());
  for (unsigned int i = 0; i < trks.size(); i++) {
    double vz = trks[i]->vz();
    if (edm::isNotFinite(vz))
      continue;
    v.add(reco::TrackBaseRef(trks[i]));
  }

  vertexes.push_back(v);

  return true;
}

bool DivisiveVertexFinder::findVertexesAlt(const reco::TrackRefVector &trks,  // input
                                           reco::VertexCollection &vertexes,
                                           const math::XYZPoint &bs) {  // output
  std::vector<PVCluster> in;
  std::pair<std::vector<PVCluster>, std::vector<const reco::Track *> > out;

  // Convert input track collection into container needed by Wolfgang's templated code
  // Need to save a map to reconvert from bare pointers, oy vey
  std::map<const reco::Track *, reco::TrackRef> mapa;
  //  std::vector< std::vector< const reco::Track* > > trkps;
  for (unsigned int i = 0; i < trks.size(); ++i) {
    double vz = trks[i]->vz();
    if (edm::isNotFinite(vz))
      continue;
    std::vector<const reco::Track *> temp;
    temp.clear();
    temp.push_back(&(*trks[i]));

    in.push_back(PVCluster(Measurement1D(trks[i]->dz(bs), trks[i]->dzError()), temp));
    mapa[temp[0]] = trks[i];
  }

  if (verbose_ > 0) {
    edm::LogInfo("DivisiveVertexFinder") << "size of input vector of clusters " << in.size();
    for (unsigned int i = 0; i < in.size(); ++i) {
      edm::LogInfo("DivisiveVertexFinder")
          << "Track " << i << " addr " << in[i].tracks()[0] << " dz " << in[i].tracks()[0]->dz(bs) << " +- "
          << in[i].tracks()[0]->dzError() << " prodID " << mapa[in[i].tracks()[0]].id() << " dz from RefTrack "
          << mapa[in[i].tracks()[0]]->dz(bs) << " +- " << mapa[in[i].tracks()[0]]->dzError();
    }
  }

  // Run the darn thing
  divmeth_.setBeamSpot(bs);
  out = divmeth_(in);

  if (verbose_ > 0)
    edm::LogInfo("DivisiveVertexFinder") << " DivisiveClusterizer1D found " << out.first.size() << " vertexes";

  // Now convert the output yet again into something we can safely store in the event
  for (unsigned int iv = 0; iv < out.first.size(); ++iv) {  // loop over output vertexes
    reco::Vertex::Error err;
    err(2, 2) = out.first[iv].position().error() * out.first[iv].position().error();

    reco::Vertex v(reco::Vertex::Point(0, 0, out.first[iv].position().value()), err, 0, 1, out.second.size());
    if (verbose_ > 0)
      edm::LogInfo("DivisiveVertexFinder")
          << " DivisiveClusterizer1D vertex " << iv << " has " << out.first[iv].tracks().size()
          << " tracks and a position of " << v.z() << " +- " << std::sqrt(v.covariance(2, 2));
    for (unsigned int itrk = 0; itrk < out.first[iv].tracks().size(); ++itrk) {
      v.add(reco::TrackBaseRef(mapa[out.first[iv].tracks()[itrk]]));
    }
    vertexes.push_back(v);  // Done with horrible conversion, save it
  }

  // Finally, sort the vertexes in decreasing sumPtSquared
  //  std::sort(vertexes.begin(), vertexes.end(), PVClusterComparer());
  std::sort(vertexes.begin(), vertexes.end(), *pvComparer_);

  return true;
}