Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:29:18

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   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   // whether to reconstruct KShorts
0061   doKShorts_ = theParameters.getParameter<bool>("doKShorts");
0062   // whether to reconstruct Lambdas
0063   doLambdas_ = theParameters.getParameter<bool>("doLambdas");
0064 
0065   // cuts on initial track selection
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   // cuts on vertex
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   // miscellaneous cuts
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   // cuts on the V0 candidate mass
0089   kShortMassCut_ = theParameters.getParameter<double>("kShortMassCut");
0090   lambdaMassCut_ = theParameters.getParameter<double>("lambdaMassCut");
0091 }
0092 
0093 // method containing the algorithm for vertex reconstruction
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   // fill vectors of TransientTracks and TrackRefs after applying preselection cuts
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   // good tracks have now been selected for vertexing
0141 
0142   // loop over tracks and vertex good charged track pairs
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         // if same-sign pairs are allowed, assign the negative and positive tracks arbitrarily
0165         negativeTrackRef = theTrackRefs[trdx1];
0166         positiveTrackRef = theTrackRefs[trdx2];
0167         negTransTkPtr = &theTransTracks[trdx1];
0168         posTransTkPtr = &theTransTracks[trdx2];
0169       }
0170 
0171       // measure distance between tracks at their closest approach
0172 
0173       //these two variables are needed to 'pin' the temporary value returned to the stack
0174       // in order to keep posState and negState from pointing to destructed objects
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       // the POCA should at least be in the sensitive volume
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       // allow for different DCA cuts depending on position of POCA
0194       if (cxPtR2 < innerOuterTkDCAThreshold_ * innerOuterTkDCAThreshold_) {
0195         if (dca > innerTkDCACut_)
0196           continue;
0197       } else {
0198         if (dca > outerTkDCACut_)
0199           continue;
0200       }
0201 
0202       // the tracks should at least point in the same quadrant
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       // calculate mPiPi
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       // Fill the vector of TransientTracks to send to KVF
0219       std::vector<reco::TransientTrack> transTracks;
0220       transTracks.reserve(2);
0221       transTracks.push_back(*posTransTkPtr);
0222       transTracks.push_back(*negTransTkPtr);
0223 
0224       // create the vertex fitter object and vertex the tracks
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       // 2D decay significance
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       // 3D decay significance
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       // make sure the vertex radius is within the inner track hit radius
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       // 2D pointing angle
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       // 3D pointing angle
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       // calculate total energy of V0 3 ways: assume it's a kShort, a Lambda, or a LambdaBar.
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       // Create momentum 4-vectors for the 3 candidate types
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       // Create the VertexCompositeCandidate object that will be stored in the Event
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       // Create daughter candidates for the VertexCompositeCandidates
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       // Store the daughter Candidates in the VertexCompositeCandidates if they pass mass cuts
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 }