Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:41

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 }