PVClusterComparer.cc

Go to the documentation of this file.

#include "RecoPixelVertexing/PixelVertexFinding/interface/PVClusterComparer.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"

#include "TMath.h"

PVClusterComparer::PVClusterComparer()
  : track_pT_min_  ( 2.5)
  , track_pT_max_  (10.)
  , track_chi2_max_(9999999.)
  , track_prob_min_(-1.)
{
  setChisquareQuantile();
}

PVClusterComparer::PVClusterComparer(double track_pt_min, double track_pt_max, double track_chi2_max, double track_prob_min)
  : track_pT_min_  (track_pt_min)
  , track_pT_max_  (track_pt_max)
  , track_chi2_max_(track_chi2_max)
  , track_prob_min_(track_prob_min)
{
  setChisquareQuantile();
}

double PVClusterComparer::pTSquaredSum(const PVCluster &v) {
  double sum=0;
  for (size_t i=0; i<v.tracks().size(); ++i) {

    double pt = v.tracks()[i]->pt();
    if (pt < track_pT_min_) continue; // Don't count tracks below track_pT_min_ (2.5 GeV)

    // RM : exclude badly reconstructed tracks from the sum
    //    if (track_prob_min_ >= 0. && track_prob_min_ <= 1.) 
    //      if (TMath::Prob(v.tracks()[i]->chi2(),v.tracks()[i]->ndof()) < track_prob_min_) continue ; 
    if (track_prob_min_ >= 0. && track_prob_min_ <= 1.) {      
      size_t ndof = v.tracks()[i]->ndof();
      if (ndof >= maxChi2_.size()) 
    updateChisquareQuantile(ndof);
      // cut on chi2 which corresponds to the configured probability
      if (v.tracks()[i]->chi2() > maxChi2_[ndof] ) continue ;   
    }
    if (v.tracks()[i]->normalizedChi2() > track_chi2_max_) continue;
    if (pt > track_pT_max_) pt = track_pT_max_;
    sum += pt*pt;
  }
  return sum;
}

double PVClusterComparer::pTSquaredSum(const reco::Vertex &v) {
  double sum=0;
  for (reco::Vertex::trackRef_iterator i=v.tracks_begin(), ie =v.tracks_end();
       i!=ie; ++i) {
    
    double pt = (*i)->pt();
    if (pt < track_pT_min_) continue; // Don't count tracks below track_pT_min_ (2.5 GeV)

    // RM : exclude badly reconstructed tracks from the sum
    //    if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
    //      if (TMath::Prob((*i)->chi2(),(*i)->ndof()) < track_prob_min_) continue ; 
    if (track_prob_min_ >= 0. && track_prob_min_ <= 1.) {      
      unsigned int ndof = (*i)->ndof();
      if (ndof >= maxChi2_.size()) 
    updateChisquareQuantile(ndof);
      // cut on chi2 which corresponds to the configured probability
      if ((*i)->chi2() > maxChi2_[ndof] ) continue ;
    }
    if ((*i)->normalizedChi2() > track_chi2_max_) continue;
    
    if (pt > track_pT_max_) pt = track_pT_max_;
    sum += pt*pt;
  }
  return sum;
}

void PVClusterComparer::setChisquareQuantile() {
  std::vector<double> maxChi2(20, 0.);
  if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
    for (size_t ndof = 0; ndof < maxChi2_.size(); ++ndof)
      // http://root.cern.ch/root/html/TMath.html#TMath:ChisquareQuantile  
      maxChi2[ndof] = TMath::ChisquareQuantile(1- track_prob_min_, ndof);
  
  maxChi2_ = maxChi2;

}

void PVClusterComparer::updateChisquareQuantile(size_t ndof) {
  
  size_t oldsize = maxChi2_.size();
  //    maxChi2_.resize(ndof+1);
  for (size_t i = oldsize; i <= ndof; ++i) {
    double chi2 = TMath::ChisquareQuantile(1- track_prob_min_, i);
    maxChi2_.push_back(chi2);
  }
}


bool PVClusterComparer::operator() (const PVCluster &v1, const PVCluster &v2) {
  return ( pTSquaredSum(v1) > pTSquaredSum(v2) );
}
bool PVClusterComparer::operator() (const reco::Vertex &v1, const reco::Vertex &v2) {
  return ( pTSquaredSum(v1) > pTSquaredSum(v2) );
}