PFTracking/src/PFDisplacedVertexHelper.cc

0001 #include "RecoParticleFlow/PFTracking/interface/PFDisplacedVertexHelper.h"
0002
0003 #include "TVector3.h"
0004 #include "DataFormats/VertexReco/interface/Vertex.h"
0005
0006 #include "TMath.h"
0007
0008 using namespace std;
0009 using namespace reco;
0010
0011 const double PFDisplacedVertexHelper::pion_mass2 = 0.0194;
0012 const double PFDisplacedVertexHelper::muon_mass2 = 0.106 * 0.106;
0013 const double PFDisplacedVertexHelper::proton_mass2 = 0.938 * 0.938;
0014
0015 //for debug only
0016 //#define PFLOW_DEBUG
0017
0018 PFDisplacedVertexHelper::PFDisplacedVertexHelper()
0019     : tracksSelector_(), vertexIdentifier_(), pvtx_(math::XYZPoint(0, 0, 0)) {}
0020
0021 PFDisplacedVertexHelper::~PFDisplacedVertexHelper() {}
0022
0023 void PFDisplacedVertexHelper::setPrimaryVertex(edm::Handle<reco::VertexCollection> mainVertexHandle,
0024                                                edm::Handle<reco::BeamSpot> beamSpotHandle) {
0025   const math::XYZPoint beamSpot = beamSpotHandle.isValid()
0026                                       ? math::XYZPoint(beamSpotHandle->x0(), beamSpotHandle->y0(), beamSpotHandle->z0())
0027                                       : math::XYZPoint(0, 0, 0);
0028
0029   // The primary vertex is taken from the refitted list,
0030   // if does not exist from the average offline beam spot position
0031   // if does not exist (0,0,0) is used
0032   pvtx_ = mainVertexHandle.isValid()
0033               ? math::XYZPoint(
0034                     mainVertexHandle->begin()->x(), mainVertexHandle->begin()->y(), mainVertexHandle->begin()->z())
0035               : beamSpot;
0036 }
0037
0038 bool PFDisplacedVertexHelper::isTrackSelected(const reco::Track& trk,
0039                                               const reco::PFDisplacedVertex::VertexTrackType vertexTrackType) const {
0040   if (!tracksSelector_.selectTracks())
0041     return true;
0042
0043   bool isGoodTrack = false;
0044
0045   bool isHighPurity = trk.quality(trk.qualityByName(tracksSelector_.quality()));
0046
0047   double nChi2 = trk.normalizedChi2();
0048   double pt = trk.pt();
0049   int nHits = trk.numberOfValidHits();
0050
0051   bool bIsPrimary = ((vertexTrackType == reco::PFDisplacedVertex::T_TO_VERTEX) ||
0052                      (vertexTrackType == reco::PFDisplacedVertex::T_MERGED));
0053
0054   if (bIsPrimary) {
0055     // Primary or merged track selection
0056     isGoodTrack = ((nChi2 > tracksSelector_.nChi2_min() && nChi2 < tracksSelector_.nChi2_max()) || isHighPurity) &&
0057                   pt > tracksSelector_.pt_min();
0058   } else {
0059     // Secondary tracks selection
0060     int nOuterHits = trk.hitPattern().numberOfLostHits(HitPattern::MISSING_OUTER_HITS);
0061
0062     double dxy = trk.dxy(pvtx_);
0063
0064     isGoodTrack = nChi2 < tracksSelector_.nChi2_max() && pt > tracksSelector_.pt_min() &&
0065                   fabs(dxy) > tracksSelector_.dxy_min() && nHits >= tracksSelector_.nHits_min() &&
0066                   nOuterHits <= tracksSelector_.nOuterHits_max();
0067   }
0068
0069   return isGoodTrack;
0070 }
0071
0072 reco::PFDisplacedVertex::VertexType PFDisplacedVertexHelper::identifyVertex(const reco::PFDisplacedVertex& v) const {
0073   if (!vertexIdentifier_.identifyVertices())
0074     return PFDisplacedVertex::ANY;
0075
0076   PFDisplacedVertex::M_Hypo massElec = PFDisplacedVertex::M_MASSLESS;
0077   PFDisplacedVertex::M_Hypo massPion = PFDisplacedVertex::M_PION;
0078
0079   math::XYZTLorentzVector mom_ee = v.secondaryMomentum(massElec, true);
0080   math::XYZTLorentzVector mom_pipi = v.secondaryMomentum(massPion, true);
0081
0082   // ===== (1) Identify fake and looper vertices ===== //
0083
0084   double ang = v.angle_io();
0085   double pt_ee = mom_ee.Pt();
0086   double eta_vtx = v.position().eta();
0087
0088   //cout << "Angle = " << ang << endl;
0089
0090   bool bDirectionFake = ang > vertexIdentifier_.angle_max();
0091   bool bLowPt = pt_ee < vertexIdentifier_.pt_min();
0092   bool bLooperEta = fabs(eta_vtx) < vertexIdentifier_.looper_eta_max();
0093
0094   bool isFake = bDirectionFake || (bLowPt && !bLooperEta);
0095   bool isLooper = !bDirectionFake && bLowPt && bLooperEta;
0096
0097   if (isFake)
0098     return PFDisplacedVertex::FAKE;
0099   if (isLooper)
0100     return PFDisplacedVertex::LOOPER;
0101
0102   // ===== (2) Identify Decays and Conversions ===== //
0103
0104   int c1 = v.originalTrack(v.refittedTracks()[0])->charge();
0105   int c2 = v.originalTrack(v.refittedTracks()[1])->charge();
0106   double mass_ee = mom_ee.M();
0107
0108   int nTracks = v.nTracks();
0109   int nSecondaryTracks = v.nSecondaryTracks();
0110   bool bPrimaryTracks = v.isTherePrimaryTracks();
0111   bool bMergedTracks = v.isThereMergedTracks();
0112
0113   bool bPair = (nTracks == nSecondaryTracks) && (nTracks == 2);
0114   bool bOpposite = (c1 * c2 < -0.1);
0115   bool bDirection = ang < vertexIdentifier_.angle_V0Conv_max();
0116   bool bConvMass = mass_ee < vertexIdentifier_.mConv_max();
0117
0118   bool bV0Conv = bPair && bOpposite && bDirection;
0119
0120   // If the basic configuration of conversions and V0 decays is respected
0121   // pair of secondary track with opposite charge and going in the right direction
0122   // the selection is then based on mass limits
0123   if (bV0Conv) {
0124     // == (2.1) Identify Conversions == //
0125
0126     bool isConv = bConvMass;
0127
0128     if (isConv)
0129       return PFDisplacedVertex::CONVERSION_LOOSE;
0130
0131     // == (2.2) Identify K0 == //
0132
0133     double mass_pipi = mom_pipi.M();
0134
0135     bool bK0Mass = mass_pipi < vertexIdentifier_.mK0_max() && mass_pipi > vertexIdentifier_.mK0_min();
0136
0137     bool isK0 = bK0Mass;
0138
0139     if (isK0)
0140       return PFDisplacedVertex::K0_DECAY;
0141
0142     // == (2.3) Identify Lambda == //
0143
0144     int lambdaKind = lambdaCP(v);
0145
0146     bool isLambda = (lambdaKind == 1);
0147     bool isLambdaBar = (lambdaKind == -1);
0148
0149     if (isLambda)
0150       return PFDisplacedVertex::LAMBDA_DECAY;
0151     if (isLambdaBar)
0152       return PFDisplacedVertex::LAMBDABAR_DECAY;
0153   }
0154
0155   // == (2.4) Identify K- and K+ ==
0156   bool bK = (nSecondaryTracks == 1) && bPrimaryTracks && !bMergedTracks && !bOpposite;
0157
0158   if (bK) {
0159     bool bKMass = isKaonMass(v);
0160
0161     bool isKPlus = bKMass && c1 > 0;
0162     bool isKMinus = bKMass && c1 < 0;
0163
0164     if (isKMinus)
0165       return PFDisplacedVertex::KMINUS_DECAY_LOOSE;
0166     if (isKPlus)
0167       return PFDisplacedVertex::KPLUS_DECAY_LOOSE;
0168   }
0169
0170   // ===== (3) Identify Nuclears, Kinks and Remaining Fakes ===== //
0171
0172   math::XYZTLorentzVector mom_prim = v.primaryMomentum(massPion, true);
0173
0174   double p_prim = mom_prim.P();
0175   double p_sec = mom_pipi.P();
0176   double pt_prim = mom_prim.Pt();
0177
0178   bool bLog = false;
0179   if (p_prim > 0.05 && p_sec > 0.05)
0180     bLog = log10(p_prim / p_sec) > vertexIdentifier_.logPrimSec_min();
0181   bool bPtMin = pt_prim > vertexIdentifier_.pt_kink_min();
0182
0183   // A vertex with at least 3 tracks is considered as high purity nuclear interaction
0184   // the only exception is K- decay into 3 prongs. To be studied.
0185   bool isNuclearHighPurity = nTracks > 2 && mass_ee > vertexIdentifier_.mNucl_min();
0186   bool isFakeHighPurity = nTracks > 2 && mass_ee < vertexIdentifier_.mNucl_min();
0187   // Two secondary tracks with some minimal tracks angular opening are still accepted
0188   // as nuclear interactions
0189   bool isNuclearLowPurity = (nTracks == nSecondaryTracks) && (nTracks == 2) && mass_ee > vertexIdentifier_.mNucl_min();
0190
0191   bool isFakeNucl = (nTracks == nSecondaryTracks) && (nTracks == 2) && mass_ee < vertexIdentifier_.mNucl_min();
0192
0193   // Kinks: 1 primary + 1 secondary is accepted only if the primary tracks
0194   // has more energy than the secondary and primary have some minimal pT
0195   // to produce a nuclear interaction
0196   bool isNuclearKink = (nSecondaryTracks == 1) && bPrimaryTracks && !bMergedTracks && bLog && bPtMin;
0197
0198   // Here some loopers may hide appearing in Particle Isolation plots. To be studied...
0199   bool isFakeKink = ((nSecondaryTracks == 1) && bMergedTracks && !bPrimaryTracks) ||
0200                     ((nSecondaryTracks == 1) && bPrimaryTracks && !bMergedTracks && (!bLog || !bPtMin));
0201
0202   if (isNuclearHighPurity)
0203     return PFDisplacedVertex::NUCL;
0204   if (isNuclearLowPurity)
0205     return PFDisplacedVertex::NUCL_LOOSE;
0206   if (isFakeKink || isFakeNucl || isFakeHighPurity)
0207     return PFDisplacedVertex::FAKE;
0208   if (isNuclearKink)
0209     return PFDisplacedVertex::NUCL_KINK;
0210
0211   return PFDisplacedVertex::ANY;
0212 }
0213
0214 int PFDisplacedVertexHelper::lambdaCP(const PFDisplacedVertex& v) const {
0215   int lambdaCP = 0;
0216
0217   vector<Track> refittedTracks = v.refittedTracks();
0218
0219   math::XYZTLorentzVector totalMomentumDcaRefit_lambda;
0220   math::XYZTLorentzVector totalMomentumDcaRefit_lambdabar;
0221
0222   reco::Track tMomentumDcaRefit_0 = refittedTracks[0];
0223   reco::Track tMomentumDcaRefit_1 = refittedTracks[1];
0224
0225   double mass2_0 = 0, mass2_1 = 0;
0226
0227   int c1 = v.originalTrack(v.refittedTracks()[1])->charge();
0228
0229   // --------------------------- lambda --------------------
0230
0231   if (c1 > 0.1)
0232     mass2_0 = pion_mass2, mass2_1 = proton_mass2;
0233   else
0234     mass2_0 = proton_mass2, mass2_1 = pion_mass2;
0235
0236   double E = sqrt(tMomentumDcaRefit_0.p() * tMomentumDcaRefit_0.p() + mass2_0);
0237
0238   math::XYZTLorentzVector momentumDcaRefit_0(
0239       tMomentumDcaRefit_0.px(), tMomentumDcaRefit_0.py(), tMomentumDcaRefit_0.pz(), E);
0240
0241   E = sqrt(tMomentumDcaRefit_1.p() * tMomentumDcaRefit_1.p() + mass2_1);
0242
0243   math::XYZTLorentzVector momentumDcaRefit_1(
0244       tMomentumDcaRefit_1.px(), tMomentumDcaRefit_1.py(), tMomentumDcaRefit_1.pz(), E);
0245
0246   totalMomentumDcaRefit_lambda = momentumDcaRefit_0 + momentumDcaRefit_1;
0247
0248   // --------------------------- anti - lambda --------------------
0249
0250   if (c1 > 0.1)
0251     mass2_1 = pion_mass2, mass2_0 = proton_mass2;
0252   else
0253     mass2_1 = proton_mass2, mass2_0 = pion_mass2;
0254
0255   E = sqrt(tMomentumDcaRefit_0.p() * tMomentumDcaRefit_0.p() + mass2_0);
0256
0257   math::XYZTLorentzVector momentumDcaRefit_01(
0258       tMomentumDcaRefit_0.px(), tMomentumDcaRefit_0.py(), tMomentumDcaRefit_0.pz(), E);
0259
0260   E = sqrt(tMomentumDcaRefit_1.p() * tMomentumDcaRefit_1.p() + mass2_1);
0261
0262   math::XYZTLorentzVector momentumDcaRefit_11(
0263       tMomentumDcaRefit_1.px(), tMomentumDcaRefit_1.py(), tMomentumDcaRefit_1.pz(), E);
0264
0265   totalMomentumDcaRefit_lambdabar = momentumDcaRefit_01 + momentumDcaRefit_11;
0266
0267   double mass_l = totalMomentumDcaRefit_lambda.M();
0268   double mass_lbar = totalMomentumDcaRefit_lambdabar.M();
0269
0270   //  cout << "mass_l = " << mass_l << " mass_lbar = " <<  mass_lbar << endl;
0271
0272   if (mass_l < mass_lbar && mass_l > vertexIdentifier_.mLambda_min() && mass_l < vertexIdentifier_.mLambda_max())
0273     lambdaCP = 1;
0274   else if (mass_lbar < mass_l && mass_lbar > vertexIdentifier_.mLambda_min() &&
0275            mass_lbar < vertexIdentifier_.mLambda_max())
0276     lambdaCP = -1;
0277   else
0278     lambdaCP = 0;
0279
0280   //cout << "lambdaCP = " << lambdaCP << endl;
0281
0282   return lambdaCP;
0283 }
0284
0285 bool PFDisplacedVertexHelper::isKaonMass(const PFDisplacedVertex& v) const {
0286   math::XYZVector trkInit = v.refittedTracks()[1].momentum();
0287   math::XYZVector trkFinal = v.refittedTracks()[0].momentum();
0288
0289   if (v.trackTypes()[0] == PFDisplacedVertex::T_TO_VERTEX) {
0290     trkInit = v.refittedTracks()[0].momentum();
0291     trkFinal = v.refittedTracks()[1].momentum();
0292   }
0293
0294   math::XYZVector trkNeutre(trkInit.x() - trkFinal.x(), trkInit.y() - trkFinal.y(), trkInit.z() - trkFinal.z());
0295
0296   double Ec = sqrt(muon_mass2 + trkFinal.Mag2());
0297   double En = sqrt(0 * 0 + trkNeutre.Mag2());
0298
0299   math::XYZTLorentzVectorD trkMuNu(trkInit.x(), trkInit.y(), trkInit.z(), Ec + En);
0300   double massMuNu = trkMuNu.M();
0301
0302   bool bKMass = massMuNu > vertexIdentifier_.mK_min() && massMuNu < vertexIdentifier_.mK_max();
0303
0304   return bKMass;
0305 }
0306
0307 void PFDisplacedVertexHelper::Dump(std::ostream& out) const {
0308   tracksSelector_.Dump();
0309   vertexIdentifier_.Dump();
0310   out << " pvtx_ = " << pvtx_ << std::endl;
0311 }