File indexing completed on 2024-04-06 12:29:18
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
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 }
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 if (useVertex_)
0053 token_vertices = iC.consumes<std::vector<reco::Vertex>>(theParameters.getParameter<edm::InputTag>("vertices"));
0054
0055 token_tracks = iC.consumes<reco::TrackCollection>(theParameters.getParameter<edm::InputTag>("trackRecoAlgorithm"));
0056 doFit_ = theParameters.getParameter<bool>("doFit");
0057 vertexFitter_ = theParameters.getParameter<bool>("vertexFitter");
0058 useRefTracks_ = theParameters.getParameter<bool>("useRefTracks");
0059
0060
0061 doKShorts_ = theParameters.getParameter<bool>("doKShorts");
0062
0063 doLambdas_ = theParameters.getParameter<bool>("doLambdas");
0064
0065
0066 tkChi2Cut_ = theParameters.getParameter<double>("tkChi2Cut");
0067 tkNHitsCut_ = theParameters.getParameter<int>("tkNHitsCut");
0068 tkPtCut_ = theParameters.getParameter<double>("tkPtCut");
0069 tkIPSigXYCut_ = theParameters.getParameter<double>("tkIPSigXYCut");
0070 tkIPSigZCut_ = theParameters.getParameter<double>("tkIPSigZCut");
0071
0072
0073 vtxChi2Cut_ = theParameters.getParameter<double>("vtxChi2Cut");
0074 vtxDecaySigXYZCut_ = theParameters.getParameter<double>("vtxDecaySigXYZCut");
0075 vtxDecaySigXYCut_ = theParameters.getParameter<double>("vtxDecaySigXYCut");
0076 vtxDecayXYCut_ = theParameters.getParameter<double>("vtxDecayXYCut");
0077 ssVtxDecayXYCut_ = theParameters.getParameter<double>("ssVtxDecayXYCut");
0078
0079 allowSS_ = theParameters.getParameter<bool>("allowSS");
0080 innerOuterTkDCAThreshold_ = theParameters.getParameter<double>("innerOuterTkDCAThreshold");
0081 innerTkDCACut_ = theParameters.getParameter<double>("innerTkDCACut");
0082 outerTkDCACut_ = theParameters.getParameter<double>("outerTkDCACut");
0083 allowWideAngleVtx_ = theParameters.getParameter<bool>("allowWideAngleVtx");
0084 mPiPiCut_ = theParameters.getParameter<double>("mPiPiCut");
0085 innerHitPosCut_ = theParameters.getParameter<double>("innerHitPosCut");
0086 cosThetaXYCut_ = theParameters.getParameter<double>("cosThetaXYCut");
0087 cosThetaXYZCut_ = theParameters.getParameter<double>("cosThetaXYZCut");
0088
0089 kShortMassCut_ = theParameters.getParameter<double>("kShortMassCut");
0090 lambdaMassCut_ = theParameters.getParameter<double>("lambdaMassCut");
0091 }
0092
0093
0094 void V0Fitter::fitAll(const edm::Event& iEvent,
0095 const edm::EventSetup& iSetup,
0096 reco::VertexCompositeCandidateCollection& theKshorts,
0097 reco::VertexCompositeCandidateCollection& theLambdas) {
0098 using std::vector;
0099
0100 edm::Handle<reco::TrackCollection> theTrackHandle;
0101 iEvent.getByToken(token_tracks, theTrackHandle);
0102 if (theTrackHandle->empty())
0103 return;
0104 const reco::TrackCollection* theTrackCollection = theTrackHandle.product();
0105
0106 edm::Handle<reco::BeamSpot> theBeamSpotHandle;
0107 iEvent.getByToken(token_beamSpot, theBeamSpotHandle);
0108 const reco::BeamSpot* theBeamSpot = theBeamSpotHandle.product();
0109 math::XYZPoint referencePos(theBeamSpot->position());
0110
0111 reco::Vertex referenceVtx;
0112 if (useVertex_) {
0113 edm::Handle<std::vector<reco::Vertex>> vertices;
0114 iEvent.getByToken(token_vertices, vertices);
0115 referenceVtx = vertices->at(0);
0116 referencePos = referenceVtx.position();
0117 }
0118
0119 const MagneticField* theMagneticField = &iSetup.getData(esTokenMF_);
0120
0121 std::vector<reco::TrackRef> theTrackRefs;
0122 std::vector<reco::TransientTrack> theTransTracks;
0123
0124
0125 for (reco::TrackCollection::const_iterator iTk = theTrackCollection->begin(); iTk != theTrackCollection->end();
0126 ++iTk) {
0127 const reco::Track* tmpTrack = &(*iTk);
0128 double ipsigXY = std::abs(tmpTrack->dxy(*theBeamSpot) / tmpTrack->dxyError());
0129 if (useVertex_)
0130 ipsigXY = std::abs(tmpTrack->dxy(referencePos) / tmpTrack->dxyError());
0131 double ipsigZ = std::abs(tmpTrack->dz(referencePos) / tmpTrack->dzError());
0132 if (tmpTrack->normalizedChi2() < tkChi2Cut_ && tmpTrack->numberOfValidHits() >= tkNHitsCut_ &&
0133 tmpTrack->pt() > tkPtCut_ && ipsigXY > tkIPSigXYCut_ && ipsigZ > tkIPSigZCut_) {
0134 reco::TrackRef tmpRef(theTrackHandle, std::distance(theTrackCollection->begin(), iTk));
0135 theTrackRefs.push_back(std::move(tmpRef));
0136 reco::TransientTrack tmpTransient(*tmpRef, theMagneticField);
0137 theTransTracks.push_back(std::move(tmpTransient));
0138 }
0139 }
0140
0141
0142
0143 for (unsigned int trdx1 = 0; trdx1 < theTrackRefs.size(); ++trdx1) {
0144 for (unsigned int trdx2 = trdx1 + 1; trdx2 < theTrackRefs.size(); ++trdx2) {
0145 reco::TrackRef positiveTrackRef;
0146 reco::TrackRef negativeTrackRef;
0147 reco::TransientTrack* posTransTkPtr = nullptr;
0148 reco::TransientTrack* negTransTkPtr = nullptr;
0149
0150 if (theTrackRefs[trdx1]->charge() < 0. && theTrackRefs[trdx2]->charge() > 0.) {
0151 negativeTrackRef = theTrackRefs[trdx1];
0152 positiveTrackRef = theTrackRefs[trdx2];
0153 negTransTkPtr = &theTransTracks[trdx1];
0154 posTransTkPtr = &theTransTracks[trdx2];
0155 } else if (theTrackRefs[trdx1]->charge() > 0. && theTrackRefs[trdx2]->charge() < 0.) {
0156 negativeTrackRef = theTrackRefs[trdx2];
0157 positiveTrackRef = theTrackRefs[trdx1];
0158 negTransTkPtr = &theTransTracks[trdx2];
0159 posTransTkPtr = &theTransTracks[trdx1];
0160 } else {
0161 if (!allowSS_)
0162 continue;
0163
0164
0165 negativeTrackRef = theTrackRefs[trdx1];
0166 positiveTrackRef = theTrackRefs[trdx2];
0167 negTransTkPtr = &theTransTracks[trdx1];
0168 posTransTkPtr = &theTransTracks[trdx2];
0169 }
0170
0171
0172
0173
0174
0175 auto const& posImpact = posTransTkPtr->impactPointTSCP();
0176 auto const& negImpact = negTransTkPtr->impactPointTSCP();
0177 if (!posImpact.isValid() || !negImpact.isValid())
0178 continue;
0179 FreeTrajectoryState const& posState = posImpact.theState();
0180 FreeTrajectoryState const& negState = negImpact.theState();
0181 ClosestApproachInRPhi cApp;
0182 cApp.calculate(posState, negState);
0183 if (!cApp.status())
0184 continue;
0185 float dca = std::abs(cApp.distance());
0186
0187
0188 GlobalPoint cxPt = cApp.crossingPoint();
0189 const double cxPtR2 = cxPt.x() * cxPt.x() + cxPt.y() * cxPt.y();
0190 if (cxPtR2 > 120. * 120. || std::abs(cxPt.z()) > 300.)
0191 continue;
0192
0193
0194 if (cxPtR2 < innerOuterTkDCAThreshold_ * innerOuterTkDCAThreshold_) {
0195 if (dca > innerTkDCACut_)
0196 continue;
0197 } else {
0198 if (dca > outerTkDCACut_)
0199 continue;
0200 }
0201
0202
0203 TrajectoryStateClosestToPoint const& posTSCP = posTransTkPtr->trajectoryStateClosestToPoint(cxPt);
0204 TrajectoryStateClosestToPoint const& negTSCP = negTransTkPtr->trajectoryStateClosestToPoint(cxPt);
0205 if (!posTSCP.isValid() || !negTSCP.isValid())
0206 continue;
0207 if (!allowWideAngleVtx_ && posTSCP.momentum().dot(negTSCP.momentum()) < 0)
0208 continue;
0209
0210
0211 double totalE = sqrt(posTSCP.momentum().mag2() + piMassSquared) + sqrt(negTSCP.momentum().mag2() + piMassSquared);
0212 double totalESq = totalE * totalE;
0213 double totalPSq = (posTSCP.momentum() + negTSCP.momentum()).mag2();
0214 double massSquared = totalESq - totalPSq;
0215 if (massSquared > mPiPiCut_ * mPiPiCut_)
0216 continue;
0217
0218
0219 std::vector<reco::TransientTrack> transTracks;
0220 transTracks.reserve(2);
0221 transTracks.push_back(*posTransTkPtr);
0222 transTracks.push_back(*negTransTkPtr);
0223
0224
0225 const GlobalError dummyError(1.0e-3, 0.0, 1.0e-3, 0.0, 0.0, 1.0e-3);
0226 TransientVertex theRecoVertex(cxPt, dummyError, transTracks, 1.0e-3);
0227 if (doFit_) {
0228 if (vertexFitter_) {
0229 KalmanVertexFitter theKalmanFitter(useRefTracks_ == 0 ? false : true);
0230 theRecoVertex = theKalmanFitter.vertex(transTracks);
0231 } else {
0232 useRefTracks_ = false;
0233 AdaptiveVertexFitter theAdaptiveFitter;
0234 theRecoVertex = theAdaptiveFitter.vertex(transTracks);
0235 }
0236 }
0237 if (!theRecoVertex.isValid())
0238 continue;
0239
0240 reco::Vertex theVtx = theRecoVertex;
0241 if (theVtx.normalizedChi2() > vtxChi2Cut_)
0242 continue;
0243 GlobalPoint vtxPos(theVtx.x(), theVtx.y(), theVtx.z());
0244
0245
0246 SMatrixSym3D totalCov = theBeamSpot->rotatedCovariance3D() + theVtx.covariance();
0247 if (useVertex_)
0248 totalCov = referenceVtx.covariance() + theVtx.covariance();
0249 SVector3 distVecXY(vtxPos.x() - referencePos.x(), vtxPos.y() - referencePos.y(), 0.);
0250 double distMagXY = ROOT::Math::Mag(distVecXY);
0251 if (distMagXY < vtxDecayXYCut_)
0252 continue;
0253 if (posTransTkPtr->charge() * negTransTkPtr->charge() > 0 && distMagXY < ssVtxDecayXYCut_)
0254 continue;
0255 double sigmaDistMagXY = sqrt(ROOT::Math::Similarity(totalCov, distVecXY)) / distMagXY;
0256 if (distMagXY < vtxDecaySigXYCut_ * sigmaDistMagXY)
0257 continue;
0258
0259
0260 if (vtxDecaySigXYZCut_ > 0.) {
0261 SVector3 distVecXYZ(
0262 vtxPos.x() - referencePos.x(), vtxPos.y() - referencePos.y(), vtxPos.z() - referencePos.z());
0263 double distMagXYZ = ROOT::Math::Mag(distVecXYZ);
0264 double sigmaDistMagXYZ = sqrt(ROOT::Math::Similarity(totalCov, distVecXYZ)) / distMagXYZ;
0265 if (distMagXYZ / sigmaDistMagXYZ < vtxDecaySigXYZCut_)
0266 continue;
0267 }
0268
0269
0270 double tkHitPosLimitSquared =
0271 (distMagXY - sigmaDistMagXY * innerHitPosCut_) * (distMagXY - sigmaDistMagXY * innerHitPosCut_);
0272 if (innerHitPosCut_ > 0. && positiveTrackRef->innerOk()) {
0273 reco::Vertex::Point posTkHitPos = positiveTrackRef->innerPosition();
0274 double posTkHitPosD2 = (posTkHitPos.x() - referencePos.x()) * (posTkHitPos.x() - referencePos.x()) +
0275 (posTkHitPos.y() - referencePos.y()) * (posTkHitPos.y() - referencePos.y());
0276 if (posTkHitPosD2 < tkHitPosLimitSquared)
0277 continue;
0278 }
0279 if (innerHitPosCut_ > 0. && negativeTrackRef->innerOk()) {
0280 reco::Vertex::Point negTkHitPos = negativeTrackRef->innerPosition();
0281 double negTkHitPosD2 = (negTkHitPos.x() - referencePos.x()) * (negTkHitPos.x() - referencePos.x()) +
0282 (negTkHitPos.y() - referencePos.y()) * (negTkHitPos.y() - referencePos.y());
0283 if (negTkHitPosD2 < tkHitPosLimitSquared)
0284 continue;
0285 }
0286
0287 std::unique_ptr<TrajectoryStateClosestToPoint> trajPlus;
0288 std::unique_ptr<TrajectoryStateClosestToPoint> trajMins;
0289 std::vector<reco::TransientTrack> theRefTracks;
0290 if (theRecoVertex.hasRefittedTracks()) {
0291 theRefTracks = theRecoVertex.refittedTracks();
0292 }
0293
0294 if (useRefTracks_ && theRefTracks.size() > 1) {
0295 reco::TransientTrack* thePositiveRefTrack = nullptr;
0296 reco::TransientTrack* theNegativeRefTrack = nullptr;
0297 for (std::vector<reco::TransientTrack>::iterator iTrack = theRefTracks.begin(); iTrack != theRefTracks.end();
0298 ++iTrack) {
0299 if (iTrack->track().charge() > 0.) {
0300 thePositiveRefTrack = &*iTrack;
0301 } else if (iTrack->track().charge() < 0.) {
0302 theNegativeRefTrack = &*iTrack;
0303 }
0304 }
0305 if (thePositiveRefTrack == nullptr || theNegativeRefTrack == nullptr)
0306 continue;
0307 trajPlus =
0308 std::make_unique<TrajectoryStateClosestToPoint>(thePositiveRefTrack->trajectoryStateClosestToPoint(vtxPos));
0309 trajMins =
0310 std::make_unique<TrajectoryStateClosestToPoint>(theNegativeRefTrack->trajectoryStateClosestToPoint(vtxPos));
0311 } else {
0312 trajPlus =
0313 std::make_unique<TrajectoryStateClosestToPoint>(posTransTkPtr->trajectoryStateClosestToPoint(vtxPos));
0314 trajMins =
0315 std::make_unique<TrajectoryStateClosestToPoint>(negTransTkPtr->trajectoryStateClosestToPoint(vtxPos));
0316 }
0317
0318 if (trajPlus.get() == nullptr || trajMins.get() == nullptr || !trajPlus->isValid() || !trajMins->isValid())
0319 continue;
0320
0321 GlobalVector positiveP(trajPlus->momentum());
0322 GlobalVector negativeP(trajMins->momentum());
0323 GlobalVector totalP(positiveP + negativeP);
0324
0325
0326 double dx = theVtx.x() - referencePos.x();
0327 double dy = theVtx.y() - referencePos.y();
0328 double px = totalP.x();
0329 double py = totalP.y();
0330 double angleXY = (dx * px + dy * py) / (sqrt(dx * dx + dy * dy) * sqrt(px * px + py * py));
0331 if (angleXY < cosThetaXYCut_)
0332 continue;
0333
0334
0335 if (cosThetaXYZCut_ > -1.) {
0336 double dz = theVtx.z() - referencePos.z();
0337 double pz = totalP.z();
0338 double angleXYZ =
0339 (dx * px + dy * py + dz * pz) / (sqrt(dx * dx + dy * dy + dz * dz) * sqrt(px * px + py * py + pz * pz));
0340 if (angleXYZ < cosThetaXYZCut_)
0341 continue;
0342 }
0343
0344
0345 double piPlusE = sqrt(positiveP.mag2() + piMassSquared);
0346 double piMinusE = sqrt(negativeP.mag2() + piMassSquared);
0347 double protonE = sqrt(positiveP.mag2() + protonMassSquared);
0348 double antiProtonE = sqrt(negativeP.mag2() + protonMassSquared);
0349 double kShortETot = piPlusE + piMinusE;
0350 double lambdaEtot = protonE + piMinusE;
0351 double lambdaBarEtot = antiProtonE + piPlusE;
0352
0353
0354 const reco::Particle::LorentzVector kShortP4(totalP.x(), totalP.y(), totalP.z(), kShortETot);
0355 const reco::Particle::LorentzVector lambdaP4(totalP.x(), totalP.y(), totalP.z(), lambdaEtot);
0356 const reco::Particle::LorentzVector lambdaBarP4(totalP.x(), totalP.y(), totalP.z(), lambdaBarEtot);
0357
0358 reco::Particle::Point vtx(theVtx.x(), theVtx.y(), theVtx.z());
0359 const reco::Vertex::CovarianceMatrix vtxCov(theVtx.covariance());
0360 double vtxChi2(theVtx.chi2());
0361 double vtxNdof(theVtx.ndof());
0362
0363
0364 reco::VertexCompositeCandidate* theKshort = nullptr;
0365 reco::VertexCompositeCandidate* theLambda = nullptr;
0366 reco::VertexCompositeCandidate* theLambdaBar = nullptr;
0367
0368 if (doKShorts_) {
0369 theKshort = new reco::VertexCompositeCandidate(0, kShortP4, vtx, vtxCov, vtxChi2, vtxNdof);
0370 }
0371 if (doLambdas_) {
0372 if (positiveP.mag2() > negativeP.mag2()) {
0373 theLambda = new reco::VertexCompositeCandidate(0, lambdaP4, vtx, vtxCov, vtxChi2, vtxNdof);
0374 } else {
0375 theLambdaBar = new reco::VertexCompositeCandidate(0, lambdaBarP4, vtx, vtxCov, vtxChi2, vtxNdof);
0376 }
0377 }
0378
0379
0380 reco::RecoChargedCandidate thePiPlusCand(
0381 1, reco::Particle::LorentzVector(positiveP.x(), positiveP.y(), positiveP.z(), piPlusE), vtx);
0382 thePiPlusCand.setTrack(positiveTrackRef);
0383
0384 reco::RecoChargedCandidate thePiMinusCand(
0385 -1, reco::Particle::LorentzVector(negativeP.x(), negativeP.y(), negativeP.z(), piMinusE), vtx);
0386 thePiMinusCand.setTrack(negativeTrackRef);
0387
0388 reco::RecoChargedCandidate theProtonCand(
0389 1, reco::Particle::LorentzVector(positiveP.x(), positiveP.y(), positiveP.z(), protonE), vtx);
0390 theProtonCand.setTrack(positiveTrackRef);
0391
0392 reco::RecoChargedCandidate theAntiProtonCand(
0393 -1, reco::Particle::LorentzVector(negativeP.x(), negativeP.y(), negativeP.z(), antiProtonE), vtx);
0394 theAntiProtonCand.setTrack(negativeTrackRef);
0395
0396 AddFourMomenta addp4;
0397
0398 if (doKShorts_) {
0399 theKshort->addDaughter(thePiPlusCand);
0400 theKshort->addDaughter(thePiMinusCand);
0401 theKshort->setPdgId(310);
0402 addp4.set(*theKshort);
0403 if (theKshort->mass() < kShortMass + kShortMassCut_ && theKshort->mass() > kShortMass - kShortMassCut_) {
0404 theKshorts.push_back(std::move(*theKshort));
0405 }
0406 }
0407 if (doLambdas_ && theLambda) {
0408 theLambda->addDaughter(theProtonCand);
0409 theLambda->addDaughter(thePiMinusCand);
0410 theLambda->setPdgId(3122);
0411 addp4.set(*theLambda);
0412 if (theLambda->mass() < lambdaMass + lambdaMassCut_ && theLambda->mass() > lambdaMass - lambdaMassCut_) {
0413 theLambdas.push_back(std::move(*theLambda));
0414 }
0415 } else if (doLambdas_ && theLambdaBar) {
0416 theLambdaBar->addDaughter(theAntiProtonCand);
0417 theLambdaBar->addDaughter(thePiPlusCand);
0418 theLambdaBar->setPdgId(-3122);
0419 addp4.set(*theLambdaBar);
0420 if (theLambdaBar->mass() < lambdaMass + lambdaMassCut_ && theLambdaBar->mass() > lambdaMass - lambdaMassCut_) {
0421 theLambdas.push_back(std::move(*theLambdaBar));
0422 }
0423 }
0424
0425 delete theKshort;
0426 delete theLambda;
0427 delete theLambdaBar;
0428 theKshort = theLambda = theLambdaBar = nullptr;
0429 }
0430 }
0431 }