PFDisplacedVertexHelper.cc

Go to the documentation of this file.

#include "RecoParticleFlow/PFTracking/interface/PFDisplacedVertexHelper.h"

#include "TVector3.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

#include "TMath.h"

using namespace std;
using namespace reco;

const double PFDisplacedVertexHelper::pion_mass2 = 0.0194;
const double PFDisplacedVertexHelper::muon_mass2 = 0.106 * 0.106;
const double PFDisplacedVertexHelper::proton_mass2 = 0.938 * 0.938;

//for debug only
//#define PFLOW_DEBUG

PFDisplacedVertexHelper::PFDisplacedVertexHelper()
    : tracksSelector_(), vertexIdentifier_(), pvtx_(math::XYZPoint(0, 0, 0)) {}

PFDisplacedVertexHelper::~PFDisplacedVertexHelper() {}

void PFDisplacedVertexHelper::setPrimaryVertex(edm::Handle<reco::VertexCollection> mainVertexHandle,
                                               edm::Handle<reco::BeamSpot> beamSpotHandle) {
  const math::XYZPoint beamSpot = beamSpotHandle.isValid()
                                      ? math::XYZPoint(beamSpotHandle->x0(), beamSpotHandle->y0(), beamSpotHandle->z0())
                                      : math::XYZPoint(0, 0, 0);

  // The primary vertex is taken from the refitted list,
  // if does not exist from the average offline beam spot position
  // if does not exist (0,0,0) is used
  pvtx_ = mainVertexHandle.isValid()
              ? math::XYZPoint(
                    mainVertexHandle->begin()->x(), mainVertexHandle->begin()->y(), mainVertexHandle->begin()->z())
              : beamSpot;
}

bool PFDisplacedVertexHelper::isTrackSelected(const reco::Track& trk,
                                              const reco::PFDisplacedVertex::VertexTrackType vertexTrackType) const {
  if (!tracksSelector_.selectTracks())
    return true;

  bool isGoodTrack = false;

  bool isHighPurity = trk.quality(trk.qualityByName(tracksSelector_.quality()));

  double nChi2 = trk.normalizedChi2();
  double pt = trk.pt();
  int nHits = trk.numberOfValidHits();

  bool bIsPrimary = ((vertexTrackType == reco::PFDisplacedVertex::T_TO_VERTEX) ||
                     (vertexTrackType == reco::PFDisplacedVertex::T_MERGED));

  if (bIsPrimary) {
    // Primary or merged track selection
    isGoodTrack = ((nChi2 > tracksSelector_.nChi2_min() && nChi2 < tracksSelector_.nChi2_max()) || isHighPurity) &&
                  pt > tracksSelector_.pt_min();
  } else {
    // Secondary tracks selection
    int nOuterHits = trk.hitPattern().numberOfLostHits(HitPattern::MISSING_OUTER_HITS);

    double dxy = trk.dxy(pvtx_);

    isGoodTrack = nChi2 < tracksSelector_.nChi2_max() && pt > tracksSelector_.pt_min() &&
                  fabs(dxy) > tracksSelector_.dxy_min() && nHits >= tracksSelector_.nHits_min() &&
                  nOuterHits <= tracksSelector_.nOuterHits_max();
  }

  return isGoodTrack;
}

reco::PFDisplacedVertex::VertexType PFDisplacedVertexHelper::identifyVertex(const reco::PFDisplacedVertex& v) const {
  if (!vertexIdentifier_.identifyVertices())
    return PFDisplacedVertex::ANY;

  PFDisplacedVertex::M_Hypo massElec = PFDisplacedVertex::M_MASSLESS;
  PFDisplacedVertex::M_Hypo massPion = PFDisplacedVertex::M_PION;

  math::XYZTLorentzVector mom_ee = v.secondaryMomentum(massElec, true);
  math::XYZTLorentzVector mom_pipi = v.secondaryMomentum(massPion, true);

  // ===== (1) Identify fake and looper vertices ===== //

  double ang = v.angle_io();
  double pt_ee = mom_ee.Pt();
  double eta_vtx = v.position().eta();

  //cout << "Angle = " << ang << endl;

  bool bDirectionFake = ang > vertexIdentifier_.angle_max();
  bool bLowPt = pt_ee < vertexIdentifier_.pt_min();
  bool bLooperEta = fabs(eta_vtx) < vertexIdentifier_.looper_eta_max();

  bool isFake = bDirectionFake || (bLowPt && !bLooperEta);
  bool isLooper = !bDirectionFake && bLowPt && bLooperEta;

  if (isFake)
    return PFDisplacedVertex::FAKE;
  if (isLooper)
    return PFDisplacedVertex::LOOPER;

  // ===== (2) Identify Decays and Conversions ===== //

  int c1 = v.originalTrack(v.refittedTracks()[0])->charge();
  int c2 = v.originalTrack(v.refittedTracks()[1])->charge();
  double mass_ee = mom_ee.M();

  int nTracks = v.nTracks();
  int nSecondaryTracks = v.nSecondaryTracks();
  bool bPrimaryTracks = v.isTherePrimaryTracks();
  bool bMergedTracks = v.isThereMergedTracks();

  bool bPair = (nTracks == nSecondaryTracks) && (nTracks == 2);
  bool bOpposite = (c1 * c2 < -0.1);
  bool bDirection = ang < vertexIdentifier_.angle_V0Conv_max();
  bool bConvMass = mass_ee < vertexIdentifier_.mConv_max();

  bool bV0Conv = bPair && bOpposite && bDirection;

  // If the basic configuration of conversions and V0 decays is respected
  // pair of secondary track with opposite charge and going in the right direction
  // the selection is then based on mass limits
  if (bV0Conv) {
    // == (2.1) Identify Conversions == //

    bool isConv = bConvMass;

    if (isConv)
      return PFDisplacedVertex::CONVERSION_LOOSE;

    // == (2.2) Identify K0 == //

    double mass_pipi = mom_pipi.M();

    bool bK0Mass = mass_pipi < vertexIdentifier_.mK0_max() && mass_pipi > vertexIdentifier_.mK0_min();

    bool isK0 = bK0Mass;

    if (isK0)
      return PFDisplacedVertex::K0_DECAY;

    // == (2.3) Identify Lambda == //

    int lambdaKind = lambdaCP(v);

    bool isLambda = (lambdaKind == 1);
    bool isLambdaBar = (lambdaKind == -1);

    if (isLambda)
      return PFDisplacedVertex::LAMBDA_DECAY;
    if (isLambdaBar)
      return PFDisplacedVertex::LAMBDABAR_DECAY;
  }

  // == (2.4) Identify K- and K+ ==
  bool bK = (nSecondaryTracks == 1) && bPrimaryTracks && !bMergedTracks && !bOpposite;

  if (bK) {
    bool bKMass = isKaonMass(v);

    bool isKPlus = bKMass && c1 > 0;
    bool isKMinus = bKMass && c1 < 0;

    if (isKMinus)
      return PFDisplacedVertex::KMINUS_DECAY_LOOSE;
    if (isKPlus)
      return PFDisplacedVertex::KPLUS_DECAY_LOOSE;
  }

  // ===== (3) Identify Nuclears, Kinks and Remaining Fakes ===== //

  math::XYZTLorentzVector mom_prim = v.primaryMomentum(massPion, true);

  double p_prim = mom_prim.P();
  double p_sec = mom_pipi.P();
  double pt_prim = mom_prim.Pt();

  bool bLog = false;
  if (p_prim > 0.05 && p_sec > 0.05)
    bLog = log10(p_prim / p_sec) > vertexIdentifier_.logPrimSec_min();
  bool bPtMin = pt_prim > vertexIdentifier_.pt_kink_min();

  // A vertex with at least 3 tracks is considered as high purity nuclear interaction
  // the only exception is K- decay into 3 prongs. To be studied.
  bool isNuclearHighPurity = nTracks > 2 && mass_ee > vertexIdentifier_.mNucl_min();
  bool isFakeHighPurity = nTracks > 2 && mass_ee < vertexIdentifier_.mNucl_min();
  // Two secondary tracks with some minimal tracks angular opening are still accepted
  // as nuclear interactions
  bool isNuclearLowPurity = (nTracks == nSecondaryTracks) && (nTracks == 2) && mass_ee > vertexIdentifier_.mNucl_min();

  bool isFakeNucl = (nTracks == nSecondaryTracks) && (nTracks == 2) && mass_ee < vertexIdentifier_.mNucl_min();

  // Kinks: 1 primary + 1 secondary is accepted only if the primary tracks
  // has more energy than the secondary and primary have some minimal pT
  // to produce a nuclear interaction
  bool isNuclearKink = (nSecondaryTracks == 1) && bPrimaryTracks && !bMergedTracks && bLog && bPtMin;

  // Here some loopers may hide appearing in Particle Isolation plots. To be studied...
  bool isFakeKink = ((nSecondaryTracks == 1) && bMergedTracks && !bPrimaryTracks) ||
                    ((nSecondaryTracks == 1) && bPrimaryTracks && !bMergedTracks && (!bLog || !bPtMin));

  if (isNuclearHighPurity)
    return PFDisplacedVertex::NUCL;
  if (isNuclearLowPurity)
    return PFDisplacedVertex::NUCL_LOOSE;
  if (isFakeKink || isFakeNucl || isFakeHighPurity)
    return PFDisplacedVertex::FAKE;
  if (isNuclearKink)
    return PFDisplacedVertex::NUCL_KINK;

  return PFDisplacedVertex::ANY;
}

int PFDisplacedVertexHelper::lambdaCP(const PFDisplacedVertex& v) const {
  int lambdaCP = 0;

  vector<Track> refittedTracks = v.refittedTracks();

  math::XYZTLorentzVector totalMomentumDcaRefit_lambda;
  math::XYZTLorentzVector totalMomentumDcaRefit_lambdabar;

  reco::Track tMomentumDcaRefit_0 = refittedTracks[0];
  reco::Track tMomentumDcaRefit_1 = refittedTracks[1];

  double mass2_0 = 0, mass2_1 = 0;

  int c1 = v.originalTrack(v.refittedTracks()[1])->charge();

  // --------------------------- lambda --------------------

  if (c1 > 0.1)
    mass2_0 = pion_mass2, mass2_1 = proton_mass2;
  else
    mass2_0 = proton_mass2, mass2_1 = pion_mass2;

  double E = sqrt(tMomentumDcaRefit_0.p() * tMomentumDcaRefit_0.p() + mass2_0);

  math::XYZTLorentzVector momentumDcaRefit_0(
      tMomentumDcaRefit_0.px(), tMomentumDcaRefit_0.py(), tMomentumDcaRefit_0.pz(), E);

  E = sqrt(tMomentumDcaRefit_1.p() * tMomentumDcaRefit_1.p() + mass2_1);

  math::XYZTLorentzVector momentumDcaRefit_1(
      tMomentumDcaRefit_1.px(), tMomentumDcaRefit_1.py(), tMomentumDcaRefit_1.pz(), E);

  totalMomentumDcaRefit_lambda = momentumDcaRefit_0 + momentumDcaRefit_1;

  // --------------------------- anti - lambda --------------------

  if (c1 > 0.1)
    mass2_1 = pion_mass2, mass2_0 = proton_mass2;
  else
    mass2_1 = proton_mass2, mass2_0 = pion_mass2;

  E = sqrt(tMomentumDcaRefit_0.p() * tMomentumDcaRefit_0.p() + mass2_0);

  math::XYZTLorentzVector momentumDcaRefit_01(
      tMomentumDcaRefit_0.px(), tMomentumDcaRefit_0.py(), tMomentumDcaRefit_0.pz(), E);

  E = sqrt(tMomentumDcaRefit_1.p() * tMomentumDcaRefit_1.p() + mass2_1);

  math::XYZTLorentzVector momentumDcaRefit_11(
      tMomentumDcaRefit_1.px(), tMomentumDcaRefit_1.py(), tMomentumDcaRefit_1.pz(), E);

  totalMomentumDcaRefit_lambdabar = momentumDcaRefit_01 + momentumDcaRefit_11;

  double mass_l = totalMomentumDcaRefit_lambda.M();
  double mass_lbar = totalMomentumDcaRefit_lambdabar.M();

  //  cout << "mass_l = " << mass_l << " mass_lbar = " <<  mass_lbar << endl;

  if (mass_l < mass_lbar && mass_l > vertexIdentifier_.mLambda_min() && mass_l < vertexIdentifier_.mLambda_max())
    lambdaCP = 1;
  else if (mass_lbar < mass_l && mass_lbar > vertexIdentifier_.mLambda_min() &&
           mass_lbar < vertexIdentifier_.mLambda_max())
    lambdaCP = -1;
  else
    lambdaCP = 0;

  //cout << "lambdaCP = " << lambdaCP << endl;

  return lambdaCP;
}

bool PFDisplacedVertexHelper::isKaonMass(const PFDisplacedVertex& v) const {
  math::XYZVector trkInit = v.refittedTracks()[1].momentum();
  math::XYZVector trkFinal = v.refittedTracks()[0].momentum();

  if (v.trackTypes()[0] == PFDisplacedVertex::T_TO_VERTEX) {
    trkInit = v.refittedTracks()[0].momentum();
    trkFinal = v.refittedTracks()[1].momentum();
  }

  math::XYZVector trkNeutre(trkInit.x() - trkFinal.x(), trkInit.y() - trkFinal.y(), trkInit.z() - trkFinal.z());

  double Ec = sqrt(muon_mass2 + trkFinal.Mag2());
  double En = sqrt(0 * 0 + trkNeutre.Mag2());

  math::XYZTLorentzVectorD trkMuNu(trkInit.x(), trkInit.y(), trkInit.z(), Ec + En);
  double massMuNu = trkMuNu.M();

  bool bKMass = massMuNu > vertexIdentifier_.mK_min() && massMuNu < vertexIdentifier_.mK_max();

  return bKMass;
}

void PFDisplacedVertexHelper::Dump(std::ostream& out) const {
  tracksSelector_.Dump();
  vertexIdentifier_.Dump();
  out << " pvtx_ = " << pvtx_ << std::endl;
}