PixelVertexFinding/src/PVClusterComparer.cc

0001 #include "RecoVertex/PixelVertexFinding/interface/PVClusterComparer.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0004
0005 #include "TMath.h"
0006
0007 PVClusterComparer::PVClusterComparer()
0008     : track_pT_min_(2.5), track_pT_max_(10.), track_chi2_max_(9999999.), track_prob_min_(-1.) {
0009   setChisquareQuantile();
0010 }
0011
0012 PVClusterComparer::PVClusterComparer(double track_pt_min,
0013                                      double track_pt_max,
0014                                      double track_chi2_max,
0015                                      double track_prob_min)
0016     : track_pT_min_(track_pt_min),
0017       track_pT_max_(track_pt_max),
0018       track_chi2_max_(track_chi2_max),
0019       track_prob_min_(track_prob_min) {
0020   setChisquareQuantile();
0021 }
0022
0023 double PVClusterComparer::pTSquaredSum(const PVCluster &v) {
0024   double sum = 0;
0025   for (size_t i = 0; i < v.tracks().size(); ++i) {
0026     double pt = v.tracks()[i]->pt();
0027     if (pt < track_pT_min_)
0028       continue;  // Don't count tracks below track_pT_min_ (2.5 GeV)
0029
0030     // RM : exclude badly reconstructed tracks from the sum
0031     //    if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
0032     //      if (TMath::Prob(v.tracks()[i]->chi2(),v.tracks()[i]->ndof()) < track_prob_min_) continue ;
0033     if (track_prob_min_ >= 0. && track_prob_min_ <= 1.) {
0034       size_t ndof = v.tracks()[i]->ndof();
0035       if (ndof >= maxChi2_.size())
0036         updateChisquareQuantile(ndof);
0037       // cut on chi2 which corresponds to the configured probability
0038       if (v.tracks()[i]->chi2() > maxChi2_[ndof])
0039         continue;
0040     }
0041     if (v.tracks()[i]->normalizedChi2() > track_chi2_max_)
0042       continue;
0043     if (pt > track_pT_max_)
0044       pt = track_pT_max_;
0045     sum += pt * pt;
0046   }
0047   return sum;
0048 }
0049
0050 double PVClusterComparer::pTSquaredSum(const reco::Vertex &v) {
0051   double sum = 0;
0052   for (reco::Vertex::trackRef_iterator i = v.tracks_begin(), ie = v.tracks_end(); i != ie; ++i) {
0053     double pt = (*i)->pt();
0054     if (pt < track_pT_min_)
0055       continue;  // Don't count tracks below track_pT_min_ (2.5 GeV)
0056
0057     // RM : exclude badly reconstructed tracks from the sum
0058     //    if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
0059     //      if (TMath::Prob((*i)->chi2(),(*i)->ndof()) < track_prob_min_) continue ;
0060     if (track_prob_min_ >= 0. && track_prob_min_ <= 1.) {
0061       unsigned int ndof = (*i)->ndof();
0062       if (ndof >= maxChi2_.size())
0063         updateChisquareQuantile(ndof);
0064       // cut on chi2 which corresponds to the configured probability
0065       if ((*i)->chi2() > maxChi2_[ndof])
0066         continue;
0067     }
0068     if ((*i)->normalizedChi2() > track_chi2_max_)
0069       continue;
0070
0071     if (pt > track_pT_max_)
0072       pt = track_pT_max_;
0073     sum += pt * pt;
0074   }
0075   return sum;
0076 }
0077
0078 void PVClusterComparer::setChisquareQuantile() {
0079   maxChi2_.clear();
0080   maxChi2_.resize(20, 0.0);
0081   if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
0082     for (size_t ndof = 0; ndof < maxChi2_.size(); ++ndof)
0083       // http://root.cern.ch/root/html/TMath.html#TMath:ChisquareQuantile
0084       maxChi2_[ndof] = TMath::ChisquareQuantile(1 - track_prob_min_, ndof);
0085 }
0086
0087 void PVClusterComparer::updateChisquareQuantile(size_t ndof) {
0088   size_t oldsize = maxChi2_.size();
0089   //    maxChi2_.resize(ndof+1);
0090   for (size_t i = oldsize; i <= ndof; ++i) {
0091     double chi2 = TMath::ChisquareQuantile(1 - track_prob_min_, i);
0092     maxChi2_.push_back(chi2);
0093   }
0094 }
0095
0096 bool PVClusterComparer::operator()(const PVCluster &v1, const PVCluster &v2) {
0097   return (pTSquaredSum(v1) > pTSquaredSum(v2));
0098 }
0099 bool PVClusterComparer::operator()(const reco::Vertex &v1, const reco::Vertex &v2) {
0100   return (pTSquaredSum(v1) > pTSquaredSum(v2));
0101 }