Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-03-26 00:00:10

0001 // -*- C++ -*-
0002 //
0003 // Package:    V0Producer
0004 // Class:      V0Fitter
0005 //
0006 /**\class V0Fitter V0Fitter.cc RecoVertex/V0Producer/src/V0Fitter.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Brian Drell
0015 //         Created:  Fri May 18 22:57:40 CEST 2007
0016 //
0017 //
0018 
0019 #include "V0Fitter.h"
0020 
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0023 #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
0024 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0025 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0026 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0027 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0028 #include <Math/Functions.h>
0029 #include <Math/SMatrix.h>
0030 #include <Math/SVector.h>
0031 #include "CommonTools/CandUtils/interface/AddFourMomenta.h"
0032 #include "DataFormats/VertexReco/interface/Vertex.h"
0033 #include <memory>
0034 #include <typeinfo>
0035 
0036 // pdg mass constants
0037 namespace {
0038   const double piMass = 0.13957018;
0039   const double piMassSquared = piMass * piMass;
0040   const double protonMass = 0.938272046;
0041   const double protonMassSquared = protonMass * protonMass;
0042   const double kShortMass = 0.497614;
0043   const double lambdaMass = 1.115683;
0044 }  // namespace
0045 
0046 typedef ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3>> SMatrixSym3D;
0047 typedef ROOT::Math::SVector<double, 3> SVector3;
0048 
0049 V0Fitter::V0Fitter(const edm::ParameterSet& theParameters, edm::ConsumesCollector&& iC) : esTokenMF_(iC.esConsumes()) {
0050   token_beamSpot = iC.consumes<reco::BeamSpot>(theParameters.getParameter<edm::InputTag>("beamSpot"));
0051   useVertex_ = theParameters.getParameter<bool>("useVertex");
0052   token_vertices = iC.consumes<std::vector<reco::Vertex>>(theParameters.getParameter<edm::InputTag>("vertices"));
0053 
0054   token_tracks = iC.consumes<reco::TrackCollection>(theParameters.getParameter<edm::InputTag>("trackRecoAlgorithm"));
0055   vertexFitter_ = theParameters.getParameter<bool>("vertexFitter");
0056   useRefTracks_ = theParameters.getParameter<bool>("useRefTracks");
0057 
0058   // whether to reconstruct KShorts
0059   doKShorts_ = theParameters.getParameter<bool>("doKShorts");
0060   // whether to reconstruct Lambdas
0061   doLambdas_ = theParameters.getParameter<bool>("doLambdas");
0062 
0063   // cuts on initial track selection
0064   tkChi2Cut_ = theParameters.getParameter<double>("tkChi2Cut");
0065   tkNHitsCut_ = theParameters.getParameter<int>("tkNHitsCut");
0066   tkPtCut_ = theParameters.getParameter<double>("tkPtCut");
0067   tkIPSigXYCut_ = theParameters.getParameter<double>("tkIPSigXYCut");
0068   tkIPSigZCut_ = theParameters.getParameter<double>("tkIPSigZCut");
0069 
0070   // cuts on vertex
0071   vtxChi2Cut_ = theParameters.getParameter<double>("vtxChi2Cut");
0072   vtxDecaySigXYZCut_ = theParameters.getParameter<double>("vtxDecaySigXYZCut");
0073   vtxDecaySigXYCut_ = theParameters.getParameter<double>("vtxDecaySigXYCut");
0074   // miscellaneous cuts
0075   tkDCACut_ = theParameters.getParameter<double>("tkDCACut");
0076   mPiPiCut_ = theParameters.getParameter<double>("mPiPiCut");
0077   innerHitPosCut_ = theParameters.getParameter<double>("innerHitPosCut");
0078   cosThetaXYCut_ = theParameters.getParameter<double>("cosThetaXYCut");
0079   cosThetaXYZCut_ = theParameters.getParameter<double>("cosThetaXYZCut");
0080   // cuts on the V0 candidate mass
0081   kShortMassCut_ = theParameters.getParameter<double>("kShortMassCut");
0082   lambdaMassCut_ = theParameters.getParameter<double>("lambdaMassCut");
0083 }
0084 
0085 // method containing the algorithm for vertex reconstruction
0086 void V0Fitter::fitAll(const edm::Event& iEvent,
0087                       const edm::EventSetup& iSetup,
0088                       reco::VertexCompositeCandidateCollection& theKshorts,
0089                       reco::VertexCompositeCandidateCollection& theLambdas) {
0090   using std::vector;
0091 
0092   edm::Handle<reco::TrackCollection> theTrackHandle;
0093   iEvent.getByToken(token_tracks, theTrackHandle);
0094   if (theTrackHandle->empty())
0095     return;
0096   const reco::TrackCollection* theTrackCollection = theTrackHandle.product();
0097 
0098   edm::Handle<reco::BeamSpot> theBeamSpotHandle;
0099   iEvent.getByToken(token_beamSpot, theBeamSpotHandle);
0100   const reco::BeamSpot* theBeamSpot = theBeamSpotHandle.product();
0101   math::XYZPoint referencePos(theBeamSpot->position());
0102 
0103   reco::Vertex referenceVtx;
0104   if (useVertex_) {
0105     edm::Handle<std::vector<reco::Vertex>> vertices;
0106     iEvent.getByToken(token_vertices, vertices);
0107     referenceVtx = vertices->at(0);
0108     referencePos = referenceVtx.position();
0109   }
0110 
0111   const MagneticField* theMagneticField = &iSetup.getData(esTokenMF_);
0112 
0113   std::vector<reco::TrackRef> theTrackRefs;
0114   std::vector<reco::TransientTrack> theTransTracks;
0115 
0116   // fill vectors of TransientTracks and TrackRefs after applying preselection cuts
0117   for (reco::TrackCollection::const_iterator iTk = theTrackCollection->begin(); iTk != theTrackCollection->end();
0118        ++iTk) {
0119     const reco::Track* tmpTrack = &(*iTk);
0120     double ipsigXY = std::abs(tmpTrack->dxy(*theBeamSpot) / tmpTrack->dxyError());
0121     if (useVertex_)
0122       ipsigXY = std::abs(tmpTrack->dxy(referencePos) / tmpTrack->dxyError());
0123     double ipsigZ = std::abs(tmpTrack->dz(referencePos) / tmpTrack->dzError());
0124     if (tmpTrack->normalizedChi2() < tkChi2Cut_ && tmpTrack->numberOfValidHits() >= tkNHitsCut_ &&
0125         tmpTrack->pt() > tkPtCut_ && ipsigXY > tkIPSigXYCut_ && ipsigZ > tkIPSigZCut_) {
0126       reco::TrackRef tmpRef(theTrackHandle, std::distance(theTrackCollection->begin(), iTk));
0127       theTrackRefs.push_back(std::move(tmpRef));
0128       reco::TransientTrack tmpTransient(*tmpRef, theMagneticField);
0129       theTransTracks.push_back(std::move(tmpTransient));
0130     }
0131   }
0132   // good tracks have now been selected for vertexing
0133 
0134   // loop over tracks and vertex good charged track pairs
0135   for (unsigned int trdx1 = 0; trdx1 < theTrackRefs.size(); ++trdx1) {
0136     for (unsigned int trdx2 = trdx1 + 1; trdx2 < theTrackRefs.size(); ++trdx2) {
0137       reco::TrackRef positiveTrackRef;
0138       reco::TrackRef negativeTrackRef;
0139       reco::TransientTrack* posTransTkPtr = nullptr;
0140       reco::TransientTrack* negTransTkPtr = nullptr;
0141 
0142       if (theTrackRefs[trdx1]->charge() < 0. && theTrackRefs[trdx2]->charge() > 0.) {
0143         negativeTrackRef = theTrackRefs[trdx1];
0144         positiveTrackRef = theTrackRefs[trdx2];
0145         negTransTkPtr = &theTransTracks[trdx1];
0146         posTransTkPtr = &theTransTracks[trdx2];
0147       } else if (theTrackRefs[trdx1]->charge() > 0. && theTrackRefs[trdx2]->charge() < 0.) {
0148         negativeTrackRef = theTrackRefs[trdx2];
0149         positiveTrackRef = theTrackRefs[trdx1];
0150         negTransTkPtr = &theTransTracks[trdx2];
0151         posTransTkPtr = &theTransTracks[trdx1];
0152       } else {
0153         continue;
0154       }
0155 
0156       // measure distance between tracks at their closest approach
0157 
0158       //these two variables are needed to 'pin' the temporary value returned to the stack
0159       // in order to keep posState and negState from pointing to destructed objects
0160       auto const& posImpact = posTransTkPtr->impactPointTSCP();
0161       auto const& negImpact = negTransTkPtr->impactPointTSCP();
0162       if (!posImpact.isValid() || !negImpact.isValid())
0163         continue;
0164       FreeTrajectoryState const& posState = posImpact.theState();
0165       FreeTrajectoryState const& negState = negImpact.theState();
0166       ClosestApproachInRPhi cApp;
0167       cApp.calculate(posState, negState);
0168       if (!cApp.status())
0169         continue;
0170       float dca = std::abs(cApp.distance());
0171       if (dca > tkDCACut_)
0172         continue;
0173 
0174       // the POCA should at least be in the sensitive volume
0175       GlobalPoint cxPt = cApp.crossingPoint();
0176       if ((cxPt.x() * cxPt.x() + cxPt.y() * cxPt.y()) > 120. * 120. || std::abs(cxPt.z()) > 300.)
0177         continue;
0178 
0179       // the tracks should at least point in the same quadrant
0180       TrajectoryStateClosestToPoint const& posTSCP = posTransTkPtr->trajectoryStateClosestToPoint(cxPt);
0181       TrajectoryStateClosestToPoint const& negTSCP = negTransTkPtr->trajectoryStateClosestToPoint(cxPt);
0182       if (!posTSCP.isValid() || !negTSCP.isValid())
0183         continue;
0184       if (posTSCP.momentum().dot(negTSCP.momentum()) < 0)
0185         continue;
0186 
0187       // calculate mPiPi
0188       double totalE = sqrt(posTSCP.momentum().mag2() + piMassSquared) + sqrt(negTSCP.momentum().mag2() + piMassSquared);
0189       double totalESq = totalE * totalE;
0190       double totalPSq = (posTSCP.momentum() + negTSCP.momentum()).mag2();
0191       double massSquared = totalESq - totalPSq;
0192       if (massSquared > mPiPiCut_ * mPiPiCut_)
0193         continue;
0194 
0195       // Fill the vector of TransientTracks to send to KVF
0196       std::vector<reco::TransientTrack> transTracks;
0197       transTracks.reserve(2);
0198       transTracks.push_back(*posTransTkPtr);
0199       transTracks.push_back(*negTransTkPtr);
0200 
0201       // create the vertex fitter object and vertex the tracks
0202       TransientVertex theRecoVertex;
0203       if (vertexFitter_) {
0204         KalmanVertexFitter theKalmanFitter(useRefTracks_ == 0 ? false : true);
0205         theRecoVertex = theKalmanFitter.vertex(transTracks);
0206       } else if (!vertexFitter_) {
0207         useRefTracks_ = false;
0208         AdaptiveVertexFitter theAdaptiveFitter;
0209         theRecoVertex = theAdaptiveFitter.vertex(transTracks);
0210       }
0211       if (!theRecoVertex.isValid())
0212         continue;
0213 
0214       reco::Vertex theVtx = theRecoVertex;
0215       if (theVtx.normalizedChi2() > vtxChi2Cut_)
0216         continue;
0217       GlobalPoint vtxPos(theVtx.x(), theVtx.y(), theVtx.z());
0218 
0219       // 2D decay significance
0220       SMatrixSym3D totalCov = theBeamSpot->rotatedCovariance3D() + theVtx.covariance();
0221       if (useVertex_)
0222         totalCov = referenceVtx.covariance() + theVtx.covariance();
0223       SVector3 distVecXY(vtxPos.x() - referencePos.x(), vtxPos.y() - referencePos.y(), 0.);
0224       double distMagXY = ROOT::Math::Mag(distVecXY);
0225       double sigmaDistMagXY = sqrt(ROOT::Math::Similarity(totalCov, distVecXY)) / distMagXY;
0226       if (distMagXY / sigmaDistMagXY < vtxDecaySigXYCut_)
0227         continue;
0228 
0229       // 3D decay significance
0230       if (vtxDecaySigXYZCut_ > 0.) {
0231         SVector3 distVecXYZ(
0232             vtxPos.x() - referencePos.x(), vtxPos.y() - referencePos.y(), vtxPos.z() - referencePos.z());
0233         double distMagXYZ = ROOT::Math::Mag(distVecXYZ);
0234         double sigmaDistMagXYZ = sqrt(ROOT::Math::Similarity(totalCov, distVecXYZ)) / distMagXYZ;
0235         if (distMagXYZ / sigmaDistMagXYZ < vtxDecaySigXYZCut_)
0236           continue;
0237       }
0238 
0239       // make sure the vertex radius is within the inner track hit radius
0240       double tkHitPosLimitSquared =
0241           (distMagXY - sigmaDistMagXY * innerHitPosCut_) * (distMagXY - sigmaDistMagXY * innerHitPosCut_);
0242       if (innerHitPosCut_ > 0. && positiveTrackRef->innerOk()) {
0243         reco::Vertex::Point posTkHitPos = positiveTrackRef->innerPosition();
0244         double posTkHitPosD2 = (posTkHitPos.x() - referencePos.x()) * (posTkHitPos.x() - referencePos.x()) +
0245                                (posTkHitPos.y() - referencePos.y()) * (posTkHitPos.y() - referencePos.y());
0246         if (posTkHitPosD2 < tkHitPosLimitSquared)
0247           continue;
0248       }
0249       if (innerHitPosCut_ > 0. && negativeTrackRef->innerOk()) {
0250         reco::Vertex::Point negTkHitPos = negativeTrackRef->innerPosition();
0251         double negTkHitPosD2 = (negTkHitPos.x() - referencePos.x()) * (negTkHitPos.x() - referencePos.x()) +
0252                                (negTkHitPos.y() - referencePos.y()) * (negTkHitPos.y() - referencePos.y());
0253         if (negTkHitPosD2 < tkHitPosLimitSquared)
0254           continue;
0255       }
0256 
0257       std::unique_ptr<TrajectoryStateClosestToPoint> trajPlus;
0258       std::unique_ptr<TrajectoryStateClosestToPoint> trajMins;
0259       std::vector<reco::TransientTrack> theRefTracks;
0260       if (theRecoVertex.hasRefittedTracks()) {
0261         theRefTracks = theRecoVertex.refittedTracks();
0262       }
0263 
0264       if (useRefTracks_ && theRefTracks.size() > 1) {
0265         reco::TransientTrack* thePositiveRefTrack = nullptr;
0266         reco::TransientTrack* theNegativeRefTrack = nullptr;
0267         for (std::vector<reco::TransientTrack>::iterator iTrack = theRefTracks.begin(); iTrack != theRefTracks.end();
0268              ++iTrack) {
0269           if (iTrack->track().charge() > 0.) {
0270             thePositiveRefTrack = &*iTrack;
0271           } else if (iTrack->track().charge() < 0.) {
0272             theNegativeRefTrack = &*iTrack;
0273           }
0274         }
0275         if (thePositiveRefTrack == nullptr || theNegativeRefTrack == nullptr)
0276           continue;
0277         trajPlus =
0278             std::make_unique<TrajectoryStateClosestToPoint>(thePositiveRefTrack->trajectoryStateClosestToPoint(vtxPos));
0279         trajMins =
0280             std::make_unique<TrajectoryStateClosestToPoint>(theNegativeRefTrack->trajectoryStateClosestToPoint(vtxPos));
0281       } else {
0282         trajPlus =
0283             std::make_unique<TrajectoryStateClosestToPoint>(posTransTkPtr->trajectoryStateClosestToPoint(vtxPos));
0284         trajMins =
0285             std::make_unique<TrajectoryStateClosestToPoint>(negTransTkPtr->trajectoryStateClosestToPoint(vtxPos));
0286       }
0287 
0288       if (trajPlus.get() == nullptr || trajMins.get() == nullptr || !trajPlus->isValid() || !trajMins->isValid())
0289         continue;
0290 
0291       GlobalVector positiveP(trajPlus->momentum());
0292       GlobalVector negativeP(trajMins->momentum());
0293       GlobalVector totalP(positiveP + negativeP);
0294 
0295       // 2D pointing angle
0296       double dx = theVtx.x() - referencePos.x();
0297       double dy = theVtx.y() - referencePos.y();
0298       double px = totalP.x();
0299       double py = totalP.y();
0300       double angleXY = (dx * px + dy * py) / (sqrt(dx * dx + dy * dy) * sqrt(px * px + py * py));
0301       if (angleXY < cosThetaXYCut_)
0302         continue;
0303 
0304       // 3D pointing angle
0305       if (cosThetaXYZCut_ > -1.) {
0306         double dz = theVtx.z() - referencePos.z();
0307         double pz = totalP.z();
0308         double angleXYZ =
0309             (dx * px + dy * py + dz * pz) / (sqrt(dx * dx + dy * dy + dz * dz) * sqrt(px * px + py * py + pz * pz));
0310         if (angleXYZ < cosThetaXYZCut_)
0311           continue;
0312       }
0313 
0314       // calculate total energy of V0 3 ways: assume it's a kShort, a Lambda, or a LambdaBar.
0315       double piPlusE = sqrt(positiveP.mag2() + piMassSquared);
0316       double piMinusE = sqrt(negativeP.mag2() + piMassSquared);
0317       double protonE = sqrt(positiveP.mag2() + protonMassSquared);
0318       double antiProtonE = sqrt(negativeP.mag2() + protonMassSquared);
0319       double kShortETot = piPlusE + piMinusE;
0320       double lambdaEtot = protonE + piMinusE;
0321       double lambdaBarEtot = antiProtonE + piPlusE;
0322 
0323       // Create momentum 4-vectors for the 3 candidate types
0324       const reco::Particle::LorentzVector kShortP4(totalP.x(), totalP.y(), totalP.z(), kShortETot);
0325       const reco::Particle::LorentzVector lambdaP4(totalP.x(), totalP.y(), totalP.z(), lambdaEtot);
0326       const reco::Particle::LorentzVector lambdaBarP4(totalP.x(), totalP.y(), totalP.z(), lambdaBarEtot);
0327 
0328       reco::Particle::Point vtx(theVtx.x(), theVtx.y(), theVtx.z());
0329       const reco::Vertex::CovarianceMatrix vtxCov(theVtx.covariance());
0330       double vtxChi2(theVtx.chi2());
0331       double vtxNdof(theVtx.ndof());
0332 
0333       // Create the VertexCompositeCandidate object that will be stored in the Event
0334       reco::VertexCompositeCandidate* theKshort = nullptr;
0335       reco::VertexCompositeCandidate* theLambda = nullptr;
0336       reco::VertexCompositeCandidate* theLambdaBar = nullptr;
0337 
0338       if (doKShorts_) {
0339         theKshort = new reco::VertexCompositeCandidate(0, kShortP4, vtx, vtxCov, vtxChi2, vtxNdof);
0340       }
0341       if (doLambdas_) {
0342         if (positiveP.mag2() > negativeP.mag2()) {
0343           theLambda = new reco::VertexCompositeCandidate(0, lambdaP4, vtx, vtxCov, vtxChi2, vtxNdof);
0344         } else {
0345           theLambdaBar = new reco::VertexCompositeCandidate(0, lambdaBarP4, vtx, vtxCov, vtxChi2, vtxNdof);
0346         }
0347       }
0348 
0349       // Create daughter candidates for the VertexCompositeCandidates
0350       reco::RecoChargedCandidate thePiPlusCand(
0351           1, reco::Particle::LorentzVector(positiveP.x(), positiveP.y(), positiveP.z(), piPlusE), vtx);
0352       thePiPlusCand.setTrack(positiveTrackRef);
0353 
0354       reco::RecoChargedCandidate thePiMinusCand(
0355           -1, reco::Particle::LorentzVector(negativeP.x(), negativeP.y(), negativeP.z(), piMinusE), vtx);
0356       thePiMinusCand.setTrack(negativeTrackRef);
0357 
0358       reco::RecoChargedCandidate theProtonCand(
0359           1, reco::Particle::LorentzVector(positiveP.x(), positiveP.y(), positiveP.z(), protonE), vtx);
0360       theProtonCand.setTrack(positiveTrackRef);
0361 
0362       reco::RecoChargedCandidate theAntiProtonCand(
0363           -1, reco::Particle::LorentzVector(negativeP.x(), negativeP.y(), negativeP.z(), antiProtonE), vtx);
0364       theAntiProtonCand.setTrack(negativeTrackRef);
0365 
0366       AddFourMomenta addp4;
0367       // Store the daughter Candidates in the VertexCompositeCandidates if they pass mass cuts
0368       if (doKShorts_) {
0369         theKshort->addDaughter(thePiPlusCand);
0370         theKshort->addDaughter(thePiMinusCand);
0371         theKshort->setPdgId(310);
0372         addp4.set(*theKshort);
0373         if (theKshort->mass() < kShortMass + kShortMassCut_ && theKshort->mass() > kShortMass - kShortMassCut_) {
0374           theKshorts.push_back(std::move(*theKshort));
0375         }
0376       }
0377       if (doLambdas_ && theLambda) {
0378         theLambda->addDaughter(theProtonCand);
0379         theLambda->addDaughter(thePiMinusCand);
0380         theLambda->setPdgId(3122);
0381         addp4.set(*theLambda);
0382         if (theLambda->mass() < lambdaMass + lambdaMassCut_ && theLambda->mass() > lambdaMass - lambdaMassCut_) {
0383           theLambdas.push_back(std::move(*theLambda));
0384         }
0385       } else if (doLambdas_ && theLambdaBar) {
0386         theLambdaBar->addDaughter(theAntiProtonCand);
0387         theLambdaBar->addDaughter(thePiPlusCand);
0388         theLambdaBar->setPdgId(-3122);
0389         addp4.set(*theLambdaBar);
0390         if (theLambdaBar->mass() < lambdaMass + lambdaMassCut_ && theLambdaBar->mass() > lambdaMass - lambdaMassCut_) {
0391           theLambdas.push_back(std::move(*theLambdaBar));
0392         }
0393       }
0394 
0395       delete theKshort;
0396       delete theLambda;
0397       delete theLambdaBar;
0398       theKshort = theLambda = theLambdaBar = nullptr;
0399     }
0400   }
0401 }