KinkFitter Class Reference

KinkFitter class to create Kink mdst's objects from reconstructed tracks. More...

#include <KinkFitter.h>

Public Member Functions
	KinkFitter (const std::string &trackFitResultsName="", const std::string &kinksName="", const std::string &recoTracksName="", const std::string &copiedRecoTracksName="RecoTracksKinkTmp")
	Constructor for the KinkFitter.
void	initializeCuts (const double vertexDistanceCut, const double vertexChi2Cut, const double precutDistance)
	Initialize the cuts which will be applied during the fit and store process.
void	setFitterMode (const unsigned char fitterMode)
	set kink fitter mode.
bool	fitAndStore (const Track *trackMother, const Track *trackDaughter, short filterFlag)
	Fit kink with cardinal hypothesis and store it if the fit was successful.

Private Member Functions
RecoTrack *	copyRecoTrackAndSplit (const RecoTrack *splitRecoTrack, const bool motherFlag, const unsigned int delta)
	Create a RecoTrack in a separate StoreArray based on one to be split.
bool	splitRecoTrack (const RecoTrack *recoTrackSplit, short &recoTrackIndexMother, short &recoTrackIndexDaughter)
	Split track into two based on |chi2/ndf - 1|.
unsigned int	combineTracksAndFit (const Track *trackMother, const Track *trackDaughter)
	Combine daughter and mother tracks in one and fit.
int	findHitPositionForReassignment (const RecoTrack *recoTrack, ROOT::Math::XYZVector &vertexPos, int direction)
	Find hit position closest to the vertex.
RecoTrack *	copyRecoTrackAndReassignCDCHits (RecoTrack *motherRecoTrack, RecoTrack *daughterRecoTrack, const bool motherFlag, const int delta)
	Copy RecoTrack to a separate StoreArray and reassign CDC hits according to delta.
bool	refitRecoTrackAfterReassign (RecoTrack *recoTrackMotherRefit, RecoTrack *recoTrackDaughterRefit, const RecoTrack *recoTrackMother, const RecoTrack *recoTrackDaughter)
	Try to fit new RecoTracks after hit reassignment.
RecoTrack *	copyRecoTrackForFlipAndRefit (const RecoTrack *recoTrack, const ROOT::Math::XYZVector &momentumSeed, const ROOT::Math::XYZVector &positionSeed, const double &timeSeed)
	Flip and refit the daughter track.
RecoTrack *	copyRecoTrackForRefit (const RecoTrack *recoTrack, const ROOT::Math::XYZVector &momentumSeed, const ROOT::Math::XYZVector &positionSeed, const double &timeSeed, const bool blockInnerStereoHits=false, const bool useAnotherFitter=false)
	Refit the daughter track blocking hits if required.
bool	isRefitImproveFilter6 (const RecoTrack *recoTrackDaughterRefit, const ROOT::Math::XYZVector &motherPosLast)
	check if the refit of filter 6 daughter tracks improves the distance between mother and daughter
bool	vertexFitWithRecoTracks (RecoTrack *recoTrackMother, RecoTrack *recoTrackDaughter, unsigned int &reassignHitStatus, ROOT::Math::XYZVector &vertexPos, double &distance, ROOT::Math::XYZVector vertexPosSeed=ROOT::Math::XYZVector(0, 0, 0))
	Fit kink vertex using RecoTrack's as inputs.
bool	extrapolateToVertex (genfit::MeasuredStateOnPlane &stMother, genfit::MeasuredStateOnPlane &stDaughter, const ROOT::Math::XYZVector &vertexPosition, unsigned int &reassignHitStatus)
	Extrapolate the fit results to the perigee to the kink vertex.
TrackFitResult *	buildTrackFitResult (RecoTrack *recoTrack, const genfit::MeasuredStateOnPlane &msop, const double Bz, const Const::ParticleType trackHypothesis)
	Build TrackFitResult of the Kink Track.
void	errMatrixForKFit (ROOT::Math::PxPyPzEVector &fourMomentum, TMatrixDSym &covMatrix6, TMatrixDSym &errMatrix7)
	Prepare the error matrix for the kFit.

Private Attributes
std::string	m_recoTracksName
	RecoTrackColName (input).
StoreArray< RecoTrack >	m_recoTracks
	RecoTrack (input)
StoreArray< TrackFitResult >	m_trackFitResults
	TrackFitResult (output).
StoreArray< Kink >	m_kinks
	Kink (output).
double	m_vertexDistanceCut
	cut on the distance at the found vertex.
double	m_vertexChi2Cut
	Chi2 cut.
double	m_precutDistance
	Preselection cut on distance between ending points of two tracks used in prefilter.
unsigned char	m_kinkFitterMode
	fitter mode from 0 to 15 written in bits:
bool	m_kinkFitterModeHitsReassignment
	fitter mode 1st bit
bool	m_kinkFitterModeFlipAndRefit
	fitter mode 2nd bit
bool	m_kinkFitterModeCombineAndFit
	fitter mode 3rd bit
bool	m_kinkFitterModeSplitTrack
	fitter mode 4th bit
StoreArray< RecoTrack >	m_copiedRecoTracks
	RecoTrack used to refit tracks.
genfit::MeasuredStateOnPlane	m_stMotherBuffer
	buffer for the MeasuredStateOnPlane of mother obtained in the vertex fit
genfit::MeasuredStateOnPlane	m_stDaughterBuffer
	buffer for the MeasuredStateOnPlane of daughter obtained in the vertex fit
RecoTrack *	m_motherKinkRecoTrackCache
	cache for the RecoTrack of mother used to find the best vertex
RecoTrack *	m_daughterKinkRecoTrackCache
	cache for the RecoTrack of daughter used to find the best vertex

Detailed Description

KinkFitter class to create Kink mdst's objects from reconstructed tracks.

To use this class, give the kinkFitter a mother and a daughter charged track and call the fitAndStore function. In case of track splitting, the mother and daughter should be the same particle.

kinkFitter.fitAndStore(B2TrackMother, B2TrackDaughter, filterFlag);

Definition at line 32 of file KinkFitter.h.

Constructor & Destructor Documentation

KinkFitter()

KinkFitter	(	const std::string &	trackFitResultsName = "",
		const std::string &	kinksName = "",
		const std::string &	recoTracksName = "",
		const std::string &	copiedRecoTracksName = "RecoTracksKinkTmp"
	)

Constructor for the KinkFitter.

Parameters

trackFitResultsName	Belle2::TrackFitResult StoreArray name.
kinksName	Kink StoreArray name.
recoTracksName	RecoTrack StoreArray name.
copiedRecoTracksName	RecoTrack StoreArray name (used for track refitting).

m_recoTracks

register m_copiedRecoTracks

relation : m_recoTracks <--> m_copiedRecoTracks

Definition at line 44 of file KinkFitter.cc.

  : m_recoTracksName(recoTracksName), m_kinkFitterMode(1)
{
  m_trackFitResults.isRequired(trackFitResultsName);
  m_kinks.registerInDataStore(kinksName, DataStore::c_WriteOut | DataStore::c_ErrorIfAlreadyRegistered);
 
  m_recoTracks.isRequired(m_recoTracksName);
 
  m_copiedRecoTracks.registerInDataStore(copiedRecoTracksName,
                                         DataStore::c_ErrorIfAlreadyRegistered);
  RecoTrack::registerRequiredRelations(m_copiedRecoTracks);
 
  m_copiedRecoTracks.registerRelationTo(m_recoTracks);
 
  B2ASSERT("Material effects not set up.  Please use SetupGenfitExtrapolationModule.",
           genfit::MaterialEffects::getInstance()->isInitialized());
  B2ASSERT("Magnetic field not set up.  Please use SetupGenfitExtrapolationModule.",
           genfit::FieldManager::getInstance()->isInitialized());
}

Member Function Documentation

buildTrackFitResult()

TrackFitResult * buildTrackFitResult ( RecoTrack *  recoTrack, const genfit::MeasuredStateOnPlane &  msop, const double  Bz, const Const::ParticleType  trackHypothesis  )

private

Build TrackFitResult of the Kink Track.

Build TrackFitResult of Kink Track.

Parameters

recoTrack	input RecoTrack
msop	MeasuredStateOnPlane extrapolated to the vertex (mother and daughter) or IP (mother)
Bz	z component of the magnetic field at vertex or IP
trackHypothesis	track hypothesis

Returns

pointer to created TrackFitResult of the Kink Track

Definition at line 115 of file KinkFitter.cc.

{
  const uint64_t hitPatternCDCInitializer = TrackBuilder::getHitPatternCDCInitializer(*recoTrack);
  const uint32_t hitPatternVXDInitializer = TrackBuilder::getHitPatternVXDInitializer(*recoTrack);
  const genfit::FitStatus* trackFitStatus = recoTrack->getTrackFitStatus();
 
  TrackFitResult* kinkTrackFitResult
    = m_trackFitResults.appendNew(ROOT::Math::XYZVector(msop.getPos()), ROOT::Math::XYZVector(msop.getMom()),
                                  msop.get6DCov(), msop.getCharge(),
                                  trackHypothesis,
                                  trackFitStatus->getPVal(),
                                  Bz, hitPatternCDCInitializer, hitPatternVXDInitializer,
                                  trackFitStatus->getNdf());
  return kinkTrackFitResult;
}

combineTracksAndFit()

unsigned int combineTracksAndFit ( const Track *  trackMother, const Track *  trackDaughter  )

private

Combine daughter and mother tracks in one and fit.

combine daughter and mother tracks and fit the result to check for clones

Parameters

trackMother	mother Track
trackDaughter	daughter Track

Returns

The filled bits of the fit result. first bit: combined pValue > mother pValue; second bit: combined pValue > daughter pValue third bit: combined NDF > daughter NDF; fourth bit: combined pValue > 10^-7; If the value > 15, the combination is assumed as failed with the following codes: 18: combined NDF < mother NDF; 19: fit failed;

Definition at line 469 of file KinkFitter.cc.

{
  RecoTrack* recoTrackMother = trackMother->getRelated<RecoTrack>(m_recoTracksName);
  RecoTrack* recoTrackDaughter = trackDaughter->getRelated<RecoTrack>(m_recoTracksName);
 
  // create a combined track by reassigning all daughter track hits to mother track
  RecoTrack* recoTrackCombinedRefit = copyRecoTrackAndReassignCDCHits(recoTrackMother,
                                      recoTrackDaughter, true, -recoTrackDaughter->getNumberOfCDCHits());
 
  // initialize fitter
  TrackFitter trackFitterDAF;
  // fit the new track
  trackFitterDAF.fit(*recoTrackCombinedRefit);
 
  // return 19 if the track fit failed
  if (!recoTrackCombinedRefit->wasFitSuccessful()) {
    B2DEBUG(29, "Refit of the combined track failed ");
    return 19;
  }
 
  // fit results of daughter, mother, and new combined tracks
  const genfit::FitStatus* motherTrackFitStatus = recoTrackMother->getTrackFitStatus();
  const genfit::FitStatus* daughterTrackFitStatus = recoTrackDaughter->getTrackFitStatus();
  const genfit::FitStatus* combinedTrackFitStatus = recoTrackCombinedRefit->getTrackFitStatus();
 
  B2DEBUG(29, "Initial mother fit result, p-value: " << motherTrackFitStatus->getPVal() << ", ndf: " <<
          motherTrackFitStatus->getNdf());
  B2DEBUG(29, "Initial daughter fit result, p-value: " << daughterTrackFitStatus->getPVal() << ", ndf: " <<
          daughterTrackFitStatus->getNdf());
  B2DEBUG(29, "Combined track fit result, p-value: " << combinedTrackFitStatus->getPVal() << ", ndf: " <<
          combinedTrackFitStatus->getNdf());
 
  // return 18 if the combined track has NDF less than mother track
  if (combinedTrackFitStatus->getNdf() < motherTrackFitStatus->getNdf())
    return 18;
 
  // filling bits according to the fit result of the combined track
  const int motherFlag = (combinedTrackFitStatus->getPVal() > motherTrackFitStatus->getPVal());
  const int daughterFlag = 2 * (combinedTrackFitStatus->getPVal() > daughterTrackFitStatus->getPVal());
  const int daughterNdfFlag = 4 * (combinedTrackFitStatus->getNdf() > daughterTrackFitStatus->getNdf());
  const int pValueFlag = 8 * (combinedTrackFitStatus->getPVal() > 0.0000001); // almost 5 sigma
 
  return motherFlag + daughterFlag + daughterNdfFlag + pValueFlag;
}

copyRecoTrackAndReassignCDCHits()

RecoTrack * copyRecoTrackAndReassignCDCHits ( RecoTrack *  motherRecoTrack, RecoTrack *  daughterRecoTrack, const bool  motherFlag, const int  delta  )

private

Copy RecoTrack to a separate StoreArray and reassign CDC hits according to delta.

Parameters

motherRecoTrack	mother RecoTrack
daughterRecoTrack	daughter RecoTrack
motherFlag	true if creating mother, false if creating daughter
delta	the position of the threshold to reassign the hits, starting from the end positive if mother hits are reassigned to daughter, negative if vice-versa

Returns

a pointer to created RecoTrack

Definition at line 193 of file KinkFitter.cc.

{
 
  // lists of CDC hits of mother and daughter tracks
  const auto motherCDCHit = motherRecoTrack->getSortedCDCHitList();
  const auto daughterCDCHit = daughterRecoTrack->getSortedCDCHitList();
 
  // initialization of helper variables
  const int deltaMother = delta > 0 ? -delta : 0;
  const int deltaDaughter = delta < 0 ? -delta : 0;
  int sortingParameterOffset = 0;
 
  // pointer to a track to be copied
  RecoTrack* recoTrackToCopy = nullptr;
  if (motherFlag) {
    recoTrackToCopy = motherRecoTrack;
  } else {
    recoTrackToCopy = daughterRecoTrack;
  }
 
  // copy recoTracks to a separate StoreArray using seed information
  RecoTrack* recoTrack = m_copiedRecoTracks.appendNew(ROOT::Math::XYZVector(recoTrackToCopy->getPositionSeed()),
                                                      ROOT::Math::XYZVector(recoTrackToCopy->getMomentumSeed()),
                                                      static_cast<short>(recoTrackToCopy->getChargeSeed()),
                                                      recoTrackToCopy->getStoreArrayNameOfPXDHits(),
                                                      recoTrackToCopy->getStoreArrayNameOfSVDHits(),
                                                      recoTrackToCopy->getStoreArrayNameOfCDCHits(),
                                                      recoTrackToCopy->getStoreArrayNameOfBKLMHits(),
                                                      recoTrackToCopy->getStoreArrayNameOfEKLMHits(),
                                                      recoTrackToCopy->getStoreArrayNameOfRecoHitInformation());
  recoTrack->setTimeSeed(recoTrackToCopy->getTimeSeed());
  recoTrack->setSeedCovariance(recoTrackToCopy->getSeedCovariance());
 
  if (motherFlag) {
 
    // copy PXD hits (we have checked in KinkFinderModule that there are no PXD hits on the other side of the track)
    for (const auto* pxdHit : recoTrackToCopy->getPXDHitList()) {
      auto recoHitInfo = recoTrackToCopy->getRecoHitInformation(pxdHit);
      recoTrack->addPXDHit(pxdHit, recoHitInfo->getSortingParameter(),
                           recoHitInfo->getFoundByTrackFinder());
    }
 
    // copy SVD hits (we have checked in KinkFinderModule that there are no SVD hits on the other side of the track)
    for (const auto* svdHit : recoTrackToCopy->getSVDHitList()) {
      auto recoHitInfo = recoTrackToCopy->getRecoHitInformation(svdHit);
      recoTrack->addSVDHit(svdHit, recoHitInfo->getSortingParameter(),
                           recoHitInfo->getFoundByTrackFinder());
    }
 
    // copy CDC hits with respect to reassignment
    for (size_t motherCDCHitIndex = 0; motherCDCHitIndex < motherCDCHit.size() + deltaMother; ++motherCDCHitIndex) {
      auto recoHitInfo = motherRecoTrack->getRecoHitInformation(motherCDCHit[motherCDCHitIndex]);
      recoTrack->addCDCHit(motherCDCHit[motherCDCHitIndex],
                           recoHitInfo->getSortingParameter(),
                           recoHitInfo->getRightLeftInformation(),
                           recoHitInfo->getFoundByTrackFinder());
 
    }
    sortingParameterOffset = recoTrack->getNumberOfTotalHits();
    for (size_t daughterCDCHitIndex = 0; daughterCDCHitIndex < static_cast<unsigned int>(deltaDaughter); ++daughterCDCHitIndex) {
      auto recoHitInfo = daughterRecoTrack->getRecoHitInformation(daughterCDCHit[daughterCDCHitIndex]);
      recoTrack->addCDCHit(daughterCDCHit[daughterCDCHitIndex],
                           recoHitInfo->getSortingParameter() + sortingParameterOffset,
                           recoHitInfo->getRightLeftInformation(),
                           recoHitInfo->getFoundByTrackFinder());
 
    }
 
    // we do not want to have KLM hits in mother track even if they exist
 
  } else {
 
    // In case of hit reassignment, we do not want to have VXD hits in the beginning of daughter track
    // even if they exist (absense of VXD hits at the end of the track was checked in KinkFinderModule)
 
    if (deltaMother)
      sortingParameterOffset = motherRecoTrack->getRecoHitInformation(motherCDCHit[motherCDCHit.size() +
                                                   deltaMother])->getSortingParameter();
    // copy CDC hits with respect to reassignment
    for (size_t motherCDCHitIndex = motherCDCHit.size() + deltaMother; motherCDCHitIndex < motherCDCHit.size(); ++motherCDCHitIndex) {
      auto recoHitInfo = motherRecoTrack->getRecoHitInformation(motherCDCHit[motherCDCHitIndex]);
      recoTrack->addCDCHit(motherCDCHit[motherCDCHitIndex],
                           recoHitInfo->getSortingParameter() - sortingParameterOffset,
                           recoHitInfo->getRightLeftInformation(),
                           recoHitInfo->getFoundByTrackFinder());
 
    }
    sortingParameterOffset = -deltaMother - deltaDaughter;
    for (size_t daughterCDCHitIndex = deltaDaughter; daughterCDCHitIndex < daughterCDCHit.size(); ++daughterCDCHitIndex) {
      auto recoHitInfo = daughterRecoTrack->getRecoHitInformation(daughterCDCHit[daughterCDCHitIndex]);
      recoTrack->addCDCHit(daughterCDCHit[daughterCDCHitIndex],
                           recoHitInfo->getSortingParameter() + sortingParameterOffset,
                           recoHitInfo->getRightLeftInformation(),
                           recoHitInfo->getFoundByTrackFinder());
 
    }
 
    // copy BKLM hits
    for (const auto* bklmHit : recoTrackToCopy->getBKLMHitList()) {
      auto recoHitInfo = recoTrackToCopy->getRecoHitInformation(bklmHit);
      recoTrack->addBKLMHit(bklmHit, recoHitInfo->getSortingParameter() + sortingParameterOffset,
                            recoHitInfo->getFoundByTrackFinder());
    }
 
    // copy EKLM hits
    for (const auto* eklmHit : recoTrackToCopy->getEKLMHitList()) {
      auto recoHitInfo = recoTrackToCopy->getRecoHitInformation(eklmHit);
      recoTrack->addEKLMHit(eklmHit, recoHitInfo->getSortingParameter() + sortingParameterOffset,
                            recoHitInfo->getFoundByTrackFinder());
    }
  }
 
  return recoTrack;
 
}

copyRecoTrackAndSplit()

RecoTrack * copyRecoTrackAndSplit ( const RecoTrack *  splitRecoTrack, const bool  motherFlag, const unsigned int  delta  )

private

Create a RecoTrack in a separate StoreArray based on one to be split.

Parameters

splitRecoTrack	RecoTrack to be split
motherFlag	true if creating mother, false if creating daughter
delta	the position of the threshold to split the hits, starting from the end

Returns

a pointer to created RecoTrack

Definition at line 515 of file KinkFitter.cc.

{
 
  // sorted list of CDC hits of track to be split
  const auto splitCDCHit = splitRecoTrack->getSortedCDCHitList();
 
  // seeds used for a copy RecoTrack
  ROOT::Math::XYZVector positionSeed(0, 0, 0);
  ROOT::Math::XYZVector momentumSeed(0, 0, 0);
  const short chargeSeed = splitRecoTrack->getTrackFitStatus()->getCharge();
  double timeSeed = 0;
  if (motherFlag) {
    positionSeed = ROOT::Math::XYZVector(splitRecoTrack->getMeasuredStateOnPlaneFromFirstHit().getPos());
    momentumSeed = ROOT::Math::XYZVector(splitRecoTrack->getMeasuredStateOnPlaneFromFirstHit().getMom());
    timeSeed = splitRecoTrack->getCardinalRepresentation()->getTime(
                 splitRecoTrack->getMeasuredStateOnPlaneFromFirstHit());
  } else {
    positionSeed = ROOT::Math::XYZVector(splitRecoTrack->getMeasuredStateOnPlaneFromLastHit().getPos());
    momentumSeed = ROOT::Math::XYZVector(splitRecoTrack->getMeasuredStateOnPlaneFromLastHit().getMom());
    timeSeed = splitRecoTrack->getCardinalRepresentation()->getTime(
                 splitRecoTrack->getMeasuredStateOnPlaneFromLastHit());
  }
 
  // copy recoTrack to a separate StoreArray using seed information
  RecoTrack* recoTrack = m_copiedRecoTracks.appendNew(positionSeed, momentumSeed, chargeSeed,
                                                      splitRecoTrack->getStoreArrayNameOfPXDHits(),
                                                      splitRecoTrack->getStoreArrayNameOfSVDHits(),
                                                      splitRecoTrack->getStoreArrayNameOfCDCHits(),
                                                      splitRecoTrack->getStoreArrayNameOfBKLMHits(),
                                                      splitRecoTrack->getStoreArrayNameOfEKLMHits(),
                                                      splitRecoTrack->getStoreArrayNameOfRecoHitInformation());
  recoTrack->setTimeSeed(timeSeed);
  recoTrack->setSeedCovariance(splitRecoTrack->getSeedCovariance());
 
  if (motherFlag) {
    // copy PXD hits (we have checked in KinkFinderModule that there no PXD hits on the other side of the track)
    for (const auto* pxdHit : splitRecoTrack->getPXDHitList()) {
      auto recoHitInfo = splitRecoTrack->getRecoHitInformation(pxdHit);
      recoTrack->addPXDHit(pxdHit, recoHitInfo->getSortingParameter(),
                           recoHitInfo->getFoundByTrackFinder());
    }
 
    // copy SVD hits (we have checked in KinkFinderModule that there no SVD hits on the other side of the track)
    for (const auto* svdHit : splitRecoTrack->getSVDHitList()) {
      auto recoHitInfo = splitRecoTrack->getRecoHitInformation(svdHit);
      recoTrack->addSVDHit(svdHit, recoHitInfo->getSortingParameter(),
                           recoHitInfo->getFoundByTrackFinder());
    }
 
 
    // copy CDC hits with respect to reassignment
    for (size_t splitCDCHitIndex = 0; splitCDCHitIndex < splitCDCHit.size() - delta; ++splitCDCHitIndex) {
      auto recoHitInfo = splitRecoTrack->getRecoHitInformation(splitCDCHit[splitCDCHitIndex]);
      recoTrack->addCDCHit(splitCDCHit[splitCDCHitIndex],
                           recoHitInfo->getSortingParameter(),
                           recoHitInfo->getRightLeftInformation(),
                           recoHitInfo->getFoundByTrackFinder());
 
    }
 
    // KLM hits are never assigned to mother track during splitting
 
  } else {
 
    // copy CDC hits with respect to reassignment
    int sortingParameterOffset = splitRecoTrack->getRecoHitInformation(splitCDCHit[splitCDCHit.size()
                                 - delta])->getSortingParameter();
    for (size_t splitCDCHitIndex = splitCDCHit.size() - delta; splitCDCHitIndex < splitCDCHit.size(); ++splitCDCHitIndex) {
      auto recoHitInfo = splitRecoTrack->getRecoHitInformation(splitCDCHit[splitCDCHitIndex]);
      recoTrack->addCDCHit(splitCDCHit[splitCDCHitIndex],
                           splitCDCHitIndex - sortingParameterOffset,
                           recoHitInfo->getRightLeftInformation(),
                           recoHitInfo->getFoundByTrackFinder());
 
    }
 
    // copy BKLM hits
    for (const auto* bklmHit : splitRecoTrack->getBKLMHitList()) {
      auto recoHitInfo = splitRecoTrack->getRecoHitInformation(bklmHit);
      recoTrack->addBKLMHit(bklmHit, recoHitInfo->getSortingParameter() - sortingParameterOffset,
                            recoHitInfo->getFoundByTrackFinder());
    }
 
    // copy EKLM hits
    for (const auto* eklmHit : splitRecoTrack->getEKLMHitList()) {
      auto recoHitInfo = splitRecoTrack->getRecoHitInformation(eklmHit);
      recoTrack->addEKLMHit(eklmHit, recoHitInfo->getSortingParameter() - sortingParameterOffset,
                            recoHitInfo->getFoundByTrackFinder());
    }
  }
 
  return recoTrack;
 
}

copyRecoTrackForFlipAndRefit()

RecoTrack * copyRecoTrackForFlipAndRefit ( const RecoTrack *  recoTrack, const ROOT::Math::XYZVector &  momentumSeed, const ROOT::Math::XYZVector &  positionSeed, const double &  timeSeed  )

private

Flip and refit the daughter track.

Parameters

recoTrack	recoTrack to flip and refit
momentumSeed	momentum seed (with the initial sign)
positionSeed	position seed
timeSeed	time seed

Returns

pointer to a new copied flipped and refitted daughter RecoTrack

Definition at line 358 of file KinkFitter.cc.

{
  // copy recoTracks to a separate StoreArray using seed information
  RecoTrack* newRecoTrack = m_copiedRecoTracks.appendNew(positionSeed, -momentumSeed,
                                                         static_cast<short>(-recoTrack->getChargeSeed()),
                                                         recoTrack->getStoreArrayNameOfPXDHits(),
                                                         recoTrack->getStoreArrayNameOfSVDHits(),
                                                         recoTrack->getStoreArrayNameOfCDCHits(),
                                                         recoTrack->getStoreArrayNameOfBKLMHits(),
                                                         recoTrack->getStoreArrayNameOfEKLMHits(),
                                                         recoTrack->getStoreArrayNameOfRecoHitInformation());
  newRecoTrack->setTimeSeed(timeSeed);
  newRecoTrack->setSeedCovariance(recoTrack->getSeedCovariance());
  newRecoTrack->addHitsFromRecoTrack(recoTrack, 0, true);
 
  // initialize fitter
  TrackFitter trackFitterKF;
  // DAF fitter usually works badly with flipped tracks, so kalmanFitter is used instead
  std::shared_ptr<genfit::KalmanFitterRefTrack> kalmanFitter = std::make_shared<genfit::KalmanFitterRefTrack>();
  kalmanFitter->setMinIterations(static_cast<unsigned int>(3));
  kalmanFitter->setMaxIterations(static_cast<unsigned int>(10));
  kalmanFitter->setMaxFailedHits(static_cast<unsigned int>(5));
  trackFitterKF.resetFitter(kalmanFitter);
  // fit the new track
  trackFitterKF.fit(*newRecoTrack);
 
  return newRecoTrack;
}

copyRecoTrackForRefit()

RecoTrack * copyRecoTrackForRefit ( const RecoTrack *  recoTrack, const ROOT::Math::XYZVector &  momentumSeed, const ROOT::Math::XYZVector &  positionSeed, const double &  timeSeed, const bool  blockInnerStereoHits = false, const bool  useAnotherFitter = false  )

private

Refit the daughter track blocking hits if required.

Parameters

recoTrack	recoTrack to refit
momentumSeed	momentum seed
positionSeed	position seed
timeSeed	time seed
blockInnerStereoHits	block the hits in the first stereo layer and all before
useAnotherFitter	use ordinary KalmanFilter

Returns

pointer to a new copied refitted daughter RecoTrack

Definition at line 391 of file KinkFitter.cc.

{
 
  // copy recoTracks to a separate StoreArray using seed information
  RecoTrack* newRecoTrack = m_copiedRecoTracks.appendNew(positionSeed, momentumSeed,
                                                         static_cast<short>(recoTrack->getChargeSeed()),
                                                         recoTrack->getStoreArrayNameOfPXDHits(),
                                                         recoTrack->getStoreArrayNameOfSVDHits(),
                                                         recoTrack->getStoreArrayNameOfCDCHits(),
                                                         recoTrack->getStoreArrayNameOfBKLMHits(),
                                                         recoTrack->getStoreArrayNameOfEKLMHits(),
                                                         recoTrack->getStoreArrayNameOfRecoHitInformation());
  newRecoTrack->setTimeSeed(timeSeed);
  newRecoTrack->setSeedCovariance(recoTrack->getSeedCovariance());
  newRecoTrack->addHitsFromRecoTrack(recoTrack);
 
  // block the hits in the first stereo layer and all before (leave at least 6 hits for fit)
  // (usually, wrong assignment of first stereo layer is responsible for wrong z coordinate)
  if (blockInnerStereoHits) {
    bool passedStereo = false;
    auto newCDCHitRefit = newRecoTrack->getSortedCDCHitList();
    for (int daughterCDCHitIndex = 0; daughterCDCHitIndex < static_cast<int>(newCDCHitRefit.size()) - 6; ++daughterCDCHitIndex) {
      if (!passedStereo && (newCDCHitRefit[daughterCDCHitIndex]->getISuperLayer() % 2 != 0))
        passedStereo = true;
      if (passedStereo && (newCDCHitRefit[daughterCDCHitIndex]->getISuperLayer() % 2 == 0))
        break;
      auto recoHitInfo = newRecoTrack->getRecoHitInformation(newCDCHitRefit[daughterCDCHitIndex]);
      recoHitInfo->setUseInFit(false);
    }
  }
 
  // fit the new track
  // initialize fitter
  TrackFitter trackFitter;
  // if useAnotherFitter true, set ordinary KalmanFilter (for filterFlag 6 it performs better than DAF)
  if (useAnotherFitter) {
    std::shared_ptr<genfit::KalmanFitterRefTrack> kalmanFitter = std::make_shared<genfit::KalmanFitterRefTrack>();
    kalmanFitter->setMinIterations(static_cast<unsigned int>(3));
    kalmanFitter->setMaxIterations(static_cast<unsigned int>(10));
    kalmanFitter->setMaxFailedHits(static_cast<unsigned int>(5));
    trackFitter.resetFitter(kalmanFitter);
 
    trackFitter.fit(*newRecoTrack);
  } else {
    trackFitter.fit(*newRecoTrack);
  }
 
  return newRecoTrack;
}

errMatrixForKFit()

void errMatrixForKFit ( ROOT::Math::PxPyPzEVector &  fourMomentum, TMatrixDSym &  covMatrix6, TMatrixDSym &  errMatrix7  )

private

Prepare the error matrix for the kFit.

Parameters

fourMomentum	four momentum of the track state to be used in kFit
covMatrix6	covariance matrix of the track state to be used in kFit
errMatrix7	error matrix of the track state to be prepared and used in kFit

Definition at line 1412 of file KinkFitter.cc.

{
 
  enum {
    c_Px, c_Py, c_Pz, c_E, c_X, c_Y, c_Z
  };
  constexpr unsigned order[] = {c_X, c_Y, c_Z, c_Px, c_Py, c_Pz};
 
  for (int i = 0; i < 6; i++) {
    for (int j = i; j < 6; j++) {
      errMatrix7(order[j], order[i]) = errMatrix7(order[i], order[j]) = covMatrix6(i, j);
    }
  }
 
  const double dEdp[] = {fourMomentum.X() / fourMomentum.E(),
                         fourMomentum.Y() / fourMomentum.E(),
                         fourMomentum.Z() / fourMomentum.E()
                        };
  constexpr unsigned componentMom[] = {c_Px, c_Py, c_Pz};
  constexpr unsigned componentPos[] = {c_X, c_Y, c_Z};
 
 
  // covariances (p,E)
  for (unsigned int comp : componentMom) {
    double covariance = 0;
    for (int k = 0; k < 3; k++) {
      covariance += errMatrix7(comp, componentMom[k]) * dEdp[k];
    }
    errMatrix7(comp, c_E) = covariance;
  }
 
  // covariances (x,E)
  for (unsigned int comp : componentPos) {
    double covariance = 0;
    for (int k = 0; k < 3; k++) {
      covariance += errMatrix7(comp, componentMom[k]) * dEdp[k];
    }
    errMatrix7(c_E, comp) = covariance;
  }
 
  // variance (E,E)
  double covariance = 0;
  for (int i = 0; i < 3; i++) {
    covariance += errMatrix7(componentMom[i], componentMom[i]) * dEdp[i] * dEdp[i];
  }
  for (int i = 0; i < 3; i++) {
    int k = (i + 1) % 3;
    covariance += 2 * errMatrix7(componentMom[i], componentMom[k]) * dEdp[i] * dEdp[k];
  }
  errMatrix7(c_E, c_E) = covariance;
 
}

extrapolateToVertex()

bool extrapolateToVertex ( genfit::MeasuredStateOnPlane &  stMother, genfit::MeasuredStateOnPlane &  stDaughter, const ROOT::Math::XYZVector &  vertexPosition, unsigned int &  reassignHitStatus  )

private

Extrapolate the fit results to the perigee to the kink vertex.

Extrapolate the states to the fitted vertex. If the vertex is inside one of the tracks, bits are set.

If the vertex is inside one of the tracks, bits in reassignHitStatus are set.

Parameters

stMother	mother MeasuredStateOnPlane
stDaughter	daughter MeasuredStateOnPlane
vertexPosition	vertex tracks to be extrapolated to
reassignHitStatus	bits to be set in case of vertex being inside one of the tracks

Returns

false if can not extrapolate one of the tracks to vertex; true otherwise

Definition at line 92 of file KinkFitter.cc.

{
  reassignHitStatus = 0;
  try {
    // extrapolate the state to the vertexPos
    // the value will be positive (negative) if the direction of the extrapolation is (counter)momentum-wise
    double extralengthMother = stMother.extrapolateToPoint(XYZToTVector(vertexPos));
    double extralengthDaughter = stDaughter.extrapolateToPoint(XYZToTVector(vertexPos));
    if (extralengthMother > 0
        && extralengthDaughter > 0) reassignHitStatus |= 0x1; // both positive means daughter hits to be reassigned to mother
    if (extralengthMother < 0
        && extralengthDaughter < 0) reassignHitStatus |= 0x2; // both negative means mother hits to be reassigned to daughter
    B2DEBUG(29, "extralengthMother=" << extralengthMother << ", extralengthDaughter=" << extralengthDaughter);
  } catch (...) {
    // Ideally, this shouldn't happen
    B2DEBUG(29, "Could not extrapolate track to vertex.");
    return false;
  }
  return true;
}

findHitPositionForReassignment()

int findHitPositionForReassignment ( const RecoTrack *  recoTrack, ROOT::Math::XYZVector &  vertexPos, int  direction  )

private

Find hit position closest to the vertex.

Parameters

recoTrack	RecoTrack
vertexPos	vertex
direction	direction (+1 for daughter track, -1 for mother track). Should be +-1

Returns

for daughter track, returns negative index of the hit, closest to the vertex for mother track, returns positive index of the hit, closest to the vertex, counting from the end of the track

Definition at line 135 of file KinkFitter.cc.

{
 
  if (direction != 1 && direction != -1) {
    B2WARNING("KinkFitter::findHitPositionForReassignment: the direction is not +-1, although should be. "
              "Set to +1 (-1) for direction > 0 (< 0).");
    if (direction > 0)
      direction = 1;
    else
      direction = -1;
  }
 
  // Helper variables to store the minimum
  double minimalDistance2 = std::numeric_limits<double>::max();
  int minimalIndex = 0;
  int riHit;
 
  // CDC Hits list to loop over
  auto cdcHits = recoTrack->getSortedCDCHitList();
 
 
  for (int cdcHitIndex = 0; cdcHitIndex < static_cast<int>(cdcHits.size()); ++cdcHitIndex) {
    if (direction > 0)
      riHit = cdcHitIndex;
    else
      riHit = static_cast<int>(cdcHits.size()) - 1 - cdcHitIndex;
 
    auto recoHitInfo = recoTrack->getRecoHitInformation(cdcHits[riHit]);
    if (!recoHitInfo->useInFit()) continue;
    try {
      const genfit::MeasuredStateOnPlane& measuredStateOnPlane = recoTrack->getMeasuredStateOnPlaneFromRecoHit(
          recoHitInfo);
      const double currentDistance2 = (ROOT::Math::XYZVector(measuredStateOnPlane.getPos()) - vertexPos).Mag2();
 
      if (currentDistance2 < minimalDistance2) {
        minimalDistance2 = currentDistance2;
        minimalIndex = cdcHitIndex;
      }
      // if it cannot find minimum in 3 iterations, stop searching
      if (cdcHitIndex - minimalIndex > 3) break;
    } catch (const NoTrackFitResult& exception) {
      B2DEBUG(29, "Can not get mSoP because of: " << exception.what());
      continue;
    } catch (const genfit::Exception& exception) {
      B2DEBUG(29, "Can not get mSoP because of: " << exception.what());
      continue;
    }
  }
  if (minimalIndex == 0)
    return -1 * direction;
  else
    return -(minimalIndex * direction);
 
}

fitAndStore()

bool fitAndStore	(	const Track *	trackMother,
		const Track *	trackDaughter,
		short	filterFlag
	)

Fit kink with cardinal hypothesis and store it if the fit was successful.

Fit and store kink.

If the corresponding flag is set, try to reassign hits between mother and daughter tracks. If the corresponding flag is set, try to flip and refit daughter track. If the corresponding flag is set, try to combine mother and daughter track and fit the resulting track. If the corresponding flag is set, try to split the track candidate selected by KinkFinderModule.

Parameters

trackMother	mother Track
trackDaughter	daughter Track (in case of splitting the same as trackMother)
filterFlag	filter with which track pair was selected Filter 1: Distance between first point of the daughter and last point of the mother < m_precutDistance (majority). Filter 2: Distance between last point of the daughter and last point of the mother < m_precutDistance (wrong daughter sign). Filter 3: Distance between the daughter Helix extrapolation to last point of the mother and last point of the mother < m_precutDistance (lost layers for daughter, second largest contribution). Filter 4: 2D distance between first point of the daughter and last point of the mother < m_precutDistance2D (bad daughter resolution recovered by hit reassignment). Filter 5: 2D distance between last point of the daughter and last point of the mother < m_precutDistance2D (bad daughter resolution and wrong daughter sign, almost no events). Filter 6: Distance between the daughter Helix extrapolation to last point of the mother and last point of the mother < m_precutDistance2D (lost layers for daughter combined with bad daughter resolution, can be recovered by refit). Filter 7: Track to split selected among mother candidates. Filter 8: Track to split selected among daughter candidates. Filter 9: Track to split selected among tracks not passing mother/daughter criteria.

Returns

true if the track pair is stored as a Kink, false in other cases

If the fitterMode requires, tries to reassign hits between daughter and mother tracks. If the fitterMode requires, tries to flip and refit the tracks.

Definition at line 828 of file KinkFitter.cc.

{
 
  // Existence of corresponding RecoTrack already checked at the module level;
  RecoTrack* recoTrackMother = trackMother->getRelated<RecoTrack>(m_recoTracksName);
  RecoTrack* recoTrackDaughter = trackDaughter->getRelated<RecoTrack>(m_recoTracksName);
 
  // Track splitting (filterFlag 7, 8, and 9)
  if (filterFlag >= 7 && filterFlag <= 9) {
    short recoTrackIndexMother = -1;
    short recoTrackIndexDaughter = -1;
    if (!splitRecoTrack(recoTrackMother, recoTrackIndexMother, recoTrackIndexDaughter))
      return false;
    recoTrackMother = m_copiedRecoTracks[recoTrackIndexMother];
    recoTrackDaughter = m_copiedRecoTracks[recoTrackIndexDaughter];
  }
 
  // Tracks selected with filterFlag from 4 to 6 are selected by 2D distance cut assuming bad z coordinate.
  // Initial refit is required for such tracks in the majority of the cases.
  // If the refit successful, use new RecoTrack for the vertex fit.
 
  // Initial refit daughter track for filterFlag 4 (mother end point and daughter start point close in 2D) and
  // filterFlag 6 (mother end point and daughter Helix extrapolation close in 2D), which do not require flipping.
  bool refitBadFlag = false;
  if (filterFlag == 4 || filterFlag == 6) {
    B2DEBUG(29, "Try to do initial refit of daughter track for filterFlag " << filterFlag);
 
    // initialize seeds for the refit
    // position of the last mother state
    ROOT::Math::XYZVector motherPosLast = ROOT::Math::XYZVector(recoTrackMother->getMeasuredStateOnPlaneFromLastHit().getPos());
    // use mother last state time as a seed for the daughter track
    double timeSeedDaughterRefit = recoTrackMother->getCardinalRepresentation()->getTime(
                                     recoTrackMother->getMeasuredStateOnPlaneFromLastHit());
    // use fitted state at the first hit for a momentum seed
    ROOT::Math::XYZVector momSeedDaughterRefit(recoTrackDaughter->getMeasuredStateOnPlaneFromFirstHit().getMom());
 
    // initialize refit conditions
    // remove hits until the first stereo layer is passed
    bool blockInnerStereoHits = false;
    if (filterFlag == 4) blockInnerStereoHits = true;
    // use ordinary KalmanFilter
    bool anotherFitter = false;
    if (filterFlag == 6) anotherFitter = true;
 
    // create a copy of the daughter track and refit it
    RecoTrack* recoTrackDaughterRefit = copyRecoTrackForRefit(recoTrackDaughter, momSeedDaughterRefit,
                                                              motherPosLast, timeSeedDaughterRefit,
                                                              blockInnerStereoHits, anotherFitter);
 
    // if the new track fit is successful, and in addition, the distance for the filterFlag 6
    // (mother end point and daughter Helix extrapolation close in 2D) is improved,
    // use it for the vertex fit
    if (recoTrackDaughterRefit->wasFitSuccessful()) {
      if (filterFlag == 4 ||
          (filterFlag == 6 && isRefitImproveFilter6(recoTrackDaughterRefit, motherPosLast)))
        recoTrackDaughter = recoTrackDaughterRefit;
      B2DEBUG(29, "Initial refit successful");
      refitBadFlag = true;
    }
  }
 
  // Flip and refit filterFlag 5 (mother end point and daughter end point close in 2D).
  if (m_kinkFitterModeFlipAndRefit && (filterFlag == 5)) {
    B2DEBUG(29, "Try to do initial flip and refit of daughter track for filterFlag " << filterFlag);
 
    // initialize seeds for the refit
    // use position of the last mother state as a seed for the daughter track
    ROOT::Math::XYZVector motherPosLast = ROOT::Math::XYZVector(recoTrackMother->getMeasuredStateOnPlaneFromLastHit().getPos());
    // use mother last state time as a seed for the daughter track
    double timeSeedDaughterFlipAndRefit = recoTrackMother->getCardinalRepresentation()->getTime(
                                            recoTrackMother->getMeasuredStateOnPlaneFromLastHit());
    // use fitted state at the last hit for a momentum seed
    ROOT::Math::XYZVector momSeedDaughterFlipAndRefit(recoTrackDaughter->getMeasuredStateOnPlaneFromLastHit().getMom());
 
    // create a copy of the daughter track, flipped and refitted
    RecoTrack* recoTrackDaughterFlipAndRefit = copyRecoTrackForFlipAndRefit(recoTrackDaughter,
                                               motherPosLast, momSeedDaughterFlipAndRefit,
                                               timeSeedDaughterFlipAndRefit);
 
    // if the new track fit is successful, use it for the vertex fit
    if (recoTrackDaughterFlipAndRefit->wasFitSuccessful()) {
      recoTrackDaughter = recoTrackDaughterFlipAndRefit;
      B2DEBUG(29, "Initial flip and refit successful");
    }
  }
 
  // fitted vertex position
  ROOT::Math::XYZVector vertexPos = ROOT::Math::XYZVector(recoTrackMother->getMeasuredStateOnPlaneFromLastHit().getPos());
 
  // flag to reassign hits, final hit to reassign, and distance at the fitted vertex
  unsigned int reassignHitStatus = 0;
  int finalHitPositionForReassignment = 0;
  double distanceAtVertex = std::numeric_limits<double>::max();
 
  // Try kink vertex fit. If the fit fails, return false immediately for all except
  // filterFlag 1 (mother end point and daughter start point close in 3D) and
  // filterFlag 3 (mother end point and daughter Helix extrapolation close in 3D).
  bool failedFitFlag = !vertexFitWithRecoTracks(recoTrackMother, recoTrackDaughter, reassignHitStatus, vertexPos, distanceAtVertex,
                                                ROOT::Math::XYZVector(recoTrackMother->getMeasuredStateOnPlaneFromLastHit().getPos()));
  if (failedFitFlag && (filterFlag != 1 && filterFlag != 3))
    return false;
 
  // If the fit fails for filterFlag 1 (mother end point and daughter start point close in 3D),
  // try to refit daughter track blocking the first stereo superlayer.
  if (failedFitFlag && filterFlag == 1) {
    B2DEBUG(29, "Try to do postVertexFit refit of daughter track for filterFlag " << filterFlag);
 
    // initialize seeds for the refit
    // position of the last mother state
    const ROOT::Math::XYZVector motherPosLast = ROOT::Math::XYZVector(recoTrackMother->getMeasuredStateOnPlaneFromLastHit().getPos());
    // use mother last state time as a seed for the daughter track
    const double timeSeedDaughterRefit = recoTrackMother->getCardinalRepresentation()->getTime(
                                           recoTrackMother->getMeasuredStateOnPlaneFromLastHit());
    // use fitted state at the first hit for a momentum seed
    const ROOT::Math::XYZVector momSeedDaughterRefit(recoTrackDaughter->getMeasuredStateOnPlaneFromFirstHit().getMom());
 
    // initialize refit conditions
    // remove hits until the first stereo layer is passed
    const bool blockInnerStereoHits = true;
    // do not use ordinary KalmanFilter
    const bool anotherFitter = false;
 
    // create a copy of the daughter track and refit it
    RecoTrack* recoTrackDaughterRefit = copyRecoTrackForRefit(recoTrackDaughter, momSeedDaughterRefit,
                                                              motherPosLast, timeSeedDaughterRefit,
                                                              blockInnerStereoHits, anotherFitter);
 
    // if the new track fit is successful, and in addition, the vertex fit is successful,
    // use a new fit and proceed. Otherwise, return false.
    if (recoTrackDaughterRefit->wasFitSuccessful() &&
        vertexFitWithRecoTracks(recoTrackMother, recoTrackDaughterRefit, reassignHitStatus, vertexPos, distanceAtVertex,
                                motherPosLast)) {
      recoTrackDaughter = recoTrackDaughterRefit;
      B2DEBUG(29, "postVertexFit refit successful");
      refitBadFlag = true;
    } else return false;
  }
 
  // If the daughter track for filterFlag 4 (mother end point and daughter start point close in 2D) or
  // filterFlag 1 (mother end point and daughter start point close in 3D) (if required) was refitted successfully,
  // reassignHitStatus may not be assigned due to blocked hits.
  // In this case, reassignment from daughter to mother is required, so we set it manually.
  if ((filterFlag == 4 || (failedFitFlag && filterFlag == 1)) && (refitBadFlag) && (reassignHitStatus == 0))
    reassignHitStatus |= 0x1;
 
  // Cache the objects used in vertex fit
  // Mother and daughter States extrapolated to the fitted vertex
  genfit::MeasuredStateOnPlane stMother = m_stMotherBuffer;
  genfit::MeasuredStateOnPlane stDaughter = m_stDaughterBuffer;
  // recoTracks used to fit the vertex
  m_motherKinkRecoTrackCache = recoTrackMother;
  m_daughterKinkRecoTrackCache = recoTrackDaughter;
 
  // if the corresponding fitterMode is used, try to reassign hits between mother and daughter tracks
  // This is done only for filterFlags, which require mother end point and daughter start point to be close to each other.
  // So filterFlags are 1, 4, and from 7 to 9 (track split).
  if (m_kinkFitterModeHitsReassignment &&
      (filterFlag == 1 || filterFlag == 4  ||
       (filterFlag >= 7 && filterFlag <= 9 && distanceAtVertex > m_vertexDistanceCut)) &&
      (reassignHitStatus != 0)) {
    B2DEBUG(29, "Start of the hits reassignment for filterFlag " << filterFlag);
 
    // initialize counter for reassigning tries
    unsigned short countReassignTries = 0;
 
    // variables to store temporary values
    int finalHitPositionForReassignmentTmp = 0;
    double distanceAtVertexTmp = std::numeric_limits<double>::max();
    ROOT::Math::XYZVector vertexPosTmp(vertexPos);
    RecoTrack* recoTrackMotherRefit = nullptr;
    RecoTrack* recoTrackDaughterRefit = nullptr;
    RecoTrack* recoTrackMotherBuffer = recoTrackMother;
    RecoTrack* recoTrackDaughterBuffer = recoTrackDaughter;
 
    // The number of tries is limited to 3
    while (reassignHitStatus != 0 && countReassignTries < 3) {
      ++countReassignTries;
      B2DEBUG(29, "Try number " << countReassignTries);
 
      // counter for failed retries
      unsigned short countBadReassignTries = 0;
 
      // find threshold hit position in the track
      // hit bellow the threshold are to be reassigned
      // positive for daughter hits reassignment to mother, negative vice-versa
      int hitPositionForReassignment = 0;
      if (reassignHitStatus & 0x1) {
 
        // daughter hits to be reassigned (negative value); direction should be 1
        hitPositionForReassignment = findHitPositionForReassignment(recoTrackDaughterBuffer, vertexPosTmp, 1);
 
        // test if the number of hits to be reassigned larger than the number of CDC hits in daughter tracks
        // minus number degree of freedom required for the fit
        if (static_cast<unsigned int>(-hitPositionForReassignment + 5) >
            recoTrackDaughterBuffer->getNumberOfCDCHits())
          break;
 
      } else if (reassignHitStatus & 0x2) {
 
        // mother hits to be reassigned (positive value); direction should be -1
        hitPositionForReassignment = findHitPositionForReassignment(recoTrackMotherBuffer, vertexPosTmp, -1);
 
        // test if the number of hits to be reassigned larger than the number of CDC hits in mother tracks
        // minus number degree of freedom required for the fit
        if (static_cast<unsigned int>(hitPositionForReassignment + 5) > recoTrackMotherBuffer->getNumberOfCDCHits()) {
 
          // if number of SVD hits is enough for the fit, reassign all the CDC hits
          if (recoTrackMotherBuffer->getNumberOfSVDHits() > 5)
            hitPositionForReassignment = recoTrackMotherBuffer->getNumberOfCDCHits();
          else
            break;
 
        }
      }
      B2DEBUG(29, "Found hit index, starting from which hits are reassigned: " << hitPositionForReassignment);
 
      // refit of the new tracks can fail when the position is too far
      // try positions closer to the end until reach it
      while (hitPositionForReassignment != 0) {
 
        // create new RecoTracks with reassigned hits in the separate StoreArray
        recoTrackMotherRefit = copyRecoTrackAndReassignCDCHits(recoTrackMotherBuffer,
                                                               recoTrackDaughterBuffer,
                                                               true, hitPositionForReassignment);
        recoTrackDaughterRefit = copyRecoTrackAndReassignCDCHits(recoTrackMotherBuffer,
                                                                 recoTrackDaughterBuffer,
                                                                 false, hitPositionForReassignment);
 
        // try to fit new RecoTracks assuming improvement of the result
        // if fit fails, try position closer to the end (no more than 5 tries)
        // if fit is successful, break the loop
        if (!refitRecoTrackAfterReassign(recoTrackMotherRefit, recoTrackDaughterRefit,
                                         recoTrackMother, recoTrackDaughter)) {
          if (hitPositionForReassignment > 0) {
            --hitPositionForReassignment;
          } else {
            ++hitPositionForReassignment;
          }
          ++countBadReassignTries;
 
          if (countBadReassignTries > 5) hitPositionForReassignment = 0;
        } else
          break;
        B2DEBUG(29, "Refit of the tracks failed, try with smaller hit index: " << hitPositionForReassignment);
      }
 
      // if hit position reaches end, the trial failed. Continue with default tracks.
      if (hitPositionForReassignment == 0) {
        B2DEBUG(29, "Reassigning of hits and refitting failed");
        break;
      }
 
      // Try to fit vertex, using previous result as a seed. If the fit fails, continue with current tracks.
      if (!vertexFitWithRecoTracks(recoTrackMotherRefit, recoTrackDaughterRefit, reassignHitStatus,
                                   vertexPosTmp, distanceAtVertexTmp, vertexPosTmp))
        break;
      recoTrackMotherBuffer = recoTrackMotherRefit;
      recoTrackDaughterBuffer = recoTrackDaughterRefit;
      finalHitPositionForReassignmentTmp += hitPositionForReassignment;
 
      // Remember the result leading to the smallest distance between tracks.
      if (distanceAtVertexTmp < m_vertexDistanceCut) {
        distanceAtVertex = distanceAtVertexTmp;
        vertexPos = vertexPosTmp;
        stMother = m_stMotherBuffer;
        stDaughter = m_stDaughterBuffer;
 
        m_motherKinkRecoTrackCache = recoTrackMotherRefit;
        m_daughterKinkRecoTrackCache = recoTrackDaughterRefit;
        finalHitPositionForReassignment = finalHitPositionForReassignmentTmp;
      }
    }
  }
 
  // If the corresponding fitterMode is used, try to flip and refit daughter track for
  // filterFlag 2 (mother end point and daughter end point close in 3D).
  if (m_kinkFitterModeFlipAndRefit && (filterFlag == 2)) {
    B2DEBUG(29, "Try to do postVertexFit flip and refit of daughter track for filterFlag " << filterFlag);
 
    // variables to store temporary values
    double distanceAtVertexTmp = std::numeric_limits<double>::max();
    ROOT::Math::XYZVector vertexPosTmp(vertexPos);
 
    // use mother last state time as a seed for the daughter track
    const double timeSeedDaughterFlipAndRefit = recoTrackMother->getCardinalRepresentation()->getTime(
                                                  recoTrackMother->getMeasuredStateOnPlaneFromLastHit());
 
    // use state at the fitted vertex for a momentum seed
    const ROOT::Math::XYZVector momSeedDaughterFlipAndRefit(stDaughter.getMom());
 
    // create a copy of the daughter track, flipped and refitted
    RecoTrack* recoTrackDaughterFlipAndRefit = copyRecoTrackForFlipAndRefit(recoTrackDaughter,
                                               vertexPos, momSeedDaughterFlipAndRefit,
                                               timeSeedDaughterFlipAndRefit);
 
    // if the vertex fit is successful and the result is improved, store it
    if (recoTrackDaughterFlipAndRefit->wasFitSuccessful() &&
        vertexFitWithRecoTracks(recoTrackMother, recoTrackDaughterFlipAndRefit, reassignHitStatus,
                                vertexPosTmp, distanceAtVertexTmp, vertexPosTmp))
      if (distanceAtVertexTmp < distanceAtVertex) {
        distanceAtVertex = distanceAtVertexTmp;
        vertexPos = vertexPosTmp;
        stMother = m_stMotherBuffer;
        stDaughter = m_stDaughterBuffer;
        m_daughterKinkRecoTrackCache = recoTrackDaughterFlipAndRefit;
 
        B2DEBUG(29, "postVertexFit flip and refit successful");
      }
  }
 
  // Try to refit daughter track for filterFlag 3 (mother end point and daughter Helix extrapolation close in 3D).
  // If it improves the vertex fit, which might even fail before, use a new result.
  if (filterFlag == 3) {
    B2DEBUG(29, "Try to do postVertexFit refit of daughter track for filterFlag " << filterFlag);
 
    // initialize seeds for the refit
    // position of the last mother state and the first daughter state
    const ROOT::Math::XYZVector motherPosLast = ROOT::Math::XYZVector(recoTrackMother->getMeasuredStateOnPlaneFromLastHit().getPos());
    // use mother last state time as a seed for the daughter track
    const double timeSeedDaughterRefit = recoTrackMother->getCardinalRepresentation()->getTime(
                                           recoTrackMother->getMeasuredStateOnPlaneFromLastHit());
    // use fitted state at the first hit for a momentum seed
    const ROOT::Math::XYZVector momSeedDaughterRefit(recoTrackDaughter->getMeasuredStateOnPlaneFromFirstHit().getMom());
 
    // variables to store temporary values
    double distanceAtVertexTmp = std::numeric_limits<double>::max();
    ROOT::Math::XYZVector vertexPosTmp(vertexPos);
 
    // initialize refit conditions
    // do not remove hits until the first stereo layer is passed
    const bool blockInnerStereoHits = false;
    // use ordinary KalmanFilter
    const bool anotherFitter = true;
 
    // create a copy of the daughter track and refit it
    RecoTrack* recoTrackDaughterRefit = copyRecoTrackForRefit(recoTrackDaughter, momSeedDaughterRefit,
                                                              motherPosLast, timeSeedDaughterRefit,
                                                              blockInnerStereoHits, anotherFitter);
 
    // if the vertex fit is successful and the result is improved, store it
    if (recoTrackDaughterRefit->wasFitSuccessful() &&
        vertexFitWithRecoTracks(recoTrackMother, recoTrackDaughterRefit, reassignHitStatus,
                                vertexPosTmp, distanceAtVertexTmp, vertexPosTmp)) {
      if ((distanceAtVertexTmp < distanceAtVertex) || failedFitFlag) {
        distanceAtVertex = distanceAtVertexTmp;
        vertexPos = vertexPosTmp;
        stMother = m_stMotherBuffer;
        stDaughter = m_stDaughterBuffer;
        m_daughterKinkRecoTrackCache = recoTrackDaughterRefit;
 
        B2DEBUG(29, "postVertexFit refit successful");
      }
    } else if (failedFitFlag) return false;
  }
 
  B2DEBUG(29, "Distance between tracks at fitted kink vertex " << distanceAtVertex);
  B2DEBUG(29, "Radius of the kink vertex " << vertexPos.Rho());
  B2DEBUG(29, "Number of reassigned hits " << finalHitPositionForReassignment);
 
  // for analysis purposes, there is no need to distinguish kinks selected with some filters,
  // so we rearrange them to simplify the output
  // value 1 is assigned to track pairs which have close endpoints
  // value 2 is assigned to track pairs which have missing layers between their endpoints,
  // so the Helix extrapolation was used
  // value 3-5 are assigned to the cases of track splitting. They define whether track to split was selected among
  // mother candidates (3), daughter candidates (4), or tracks not passing any of these two criteria (5)
  short filterFlagToStore = 0;
  switch (filterFlag) {
    case 1:
    case 2:
    case 4:
    case 5:
      filterFlagToStore = 1;
      break;
    case 3:
    case 6:
      filterFlagToStore = 2;
      break;
    case 7:
      filterFlagToStore = 3;
      break;
    case 8:
      filterFlagToStore = 4;
      break;
    case 9:
      filterFlagToStore = 5;
  }
 
  // test the distance cut and remove pairs that do not pass the criteria
  // for split tracks, we do not remove the candidates, but fill the second digit of the flag to store
  if (distanceAtVertex > m_vertexDistanceCut) {
    if (filterFlag < 7)
      return false;
    else
      filterFlagToStore += 10;
  }
 
  // check if the fitted vertex is inside CDC or just after SVD
  if (vertexPos.Rho() < 14)
    return false;
 
  // extrapolate the mother state to IP (B2Vector3D(0., 0., 0.) beam spot, and B2Vector3D(0., 0., 1.) beam axis)
  genfit::MeasuredStateOnPlane stMotherIP = recoTrackMother->getMeasuredStateOnPlaneFromFirstHit();
  try {
    stMotherIP.extrapolateToLine(B2Vector3D(0., 0., 0.), B2Vector3D(0., 0., 1.));
  } catch (...) {
    B2DEBUG(29, "Could not extrapolate mother track to IP.");
  }
 
  // magnetic field at the fitted vertex and IP
  const double BzVtx = BFieldManager::getFieldInTesla(vertexPos).Z();
  const double BzIP = BFieldManager::getFieldInTesla({0, 0, 0}).Z();
 
  // prepare TrackFitResults for mother at IP and fitted vertex and for daughter at fitted vertex
  TrackFitResult* tfrMotherIP = buildTrackFitResult(m_motherKinkRecoTrackCache, stMotherIP, BzIP, Const::pion);
  TrackFitResult* tfrMotherVtx = buildTrackFitResult(m_motherKinkRecoTrackCache, stMother, BzVtx, Const::pion);
  TrackFitResult* tfrDaughterVtx = buildTrackFitResult(m_daughterKinkRecoTrackCache, stDaughter, BzVtx, Const::pion);
 
  // Try to combine tracks and fit them to find clones (excluding split tracks).
  // The result is written in the second and third digits of filter flag.
  if (m_kinkFitterModeCombineAndFit && (filterFlag < 7)) {
    unsigned int combinedFitFlag = combineTracksAndFit(trackMother, trackDaughter);
    filterFlagToStore += combinedFitFlag * 10;
  }
 
  // write to the filter flag number of reassigned hits (minus for daughter to mother, plus vice-versa)
  // since the type of flag to store is short, we are limited by +-32768.
  // The number of reassigned hits is rarely exceeds 32, so for that cases we feel the corresponding digits
  // with maximum available value of 32
  if (abs(finalHitPositionForReassignment) < 32) {
    if (finalHitPositionForReassignment >= 0)
      filterFlagToStore += finalHitPositionForReassignment * 1000;
    else {
      filterFlagToStore *= -1;
      filterFlagToStore += finalHitPositionForReassignment * 1000;
    }
  } else {
    filterFlagToStore += 32 * 1000;
  }
 
  // save the kink to the StoreArray
  m_kinks.appendNew(std::make_pair(trackMother, std::make_pair(tfrMotherIP, tfrMotherVtx)),
                    std::make_pair(trackDaughter, tfrDaughterVtx),
                    vertexPos.X(), vertexPos.Y(), vertexPos.Z(), filterFlagToStore);
 
 
  return true;
}

initializeCuts()

void initializeCuts	(	const double	vertexDistanceCut,
		const double	vertexChi2Cut,
		const double	precutDistance
	)

Initialize the cuts which will be applied during the fit and store process.

Parameters

vertexDistanceCut	Cut on distance between tracks at the Kink vertex.
vertexChi2Cut	Cut on Chi2 for the Kink vertex.
precutDistance	Preselection cut on distance between ending points of two tracks used in KinkFinderModule.

Definition at line 82 of file KinkFitter.cc.

{
  m_vertexDistanceCut = vertexDistanceCut;
  m_vertexChi2Cut = vertexChi2Cut;
  m_precutDistance = precutDistance;
}

isRefitImproveFilter6()

bool isRefitImproveFilter6 ( const RecoTrack *  recoTrackDaughterRefit, const ROOT::Math::XYZVector &  motherPosLast  )

private

check if the refit of filter 6 daughter tracks improves the distance between mother and daughter

check if the refit of daughter track for filterFlag 6 (mother end point and daughter Helix extrapolation close in 2D) improves the 3D distance between mother and daughter

Parameters

recoTrackDaughterRefit	refitted daughter recoTrack
motherPosLast	position of the mother state at last hit

Returns

true if the refit of filter 6 daughter tracks improves the distance between mother and daughter; false otherwise

Definition at line 447 of file KinkFitter.cc.

{
  // get the values near the mother last point
  ROOT::Math::XYZVector daughterPosClosestToMotherPosLast = ROOT::Math::XYZVector(
                                                              recoTrackDaughterRefit->getMeasuredStateOnPlaneFromFirstHit().getPos());
  ROOT::Math::XYZVector daughterMomClosestToMotherPosLast = ROOT::Math::XYZVector(
                                                              recoTrackDaughterRefit->getMeasuredStateOnPlaneFromFirstHit().getMom());
  const double Bz = BFieldManager::getFieldInTesla(daughterPosClosestToMotherPosLast).Z();
  // daughter Helix with move to mother last point
  Helix daughterHelixClosestToMotherPosLast(daughterPosClosestToMotherPosLast,
                                            daughterMomClosestToMotherPosLast,
                                            static_cast<short>(recoTrackDaughterRefit->getTrackFitStatus()->getCharge()),
                                            Bz);
  daughterHelixClosestToMotherPosLast.passiveMoveBy(motherPosLast);
 
  // check if the 3D distance passes loose criteria (in default track fit used in KinkFinderModule this test is failed)
  return ((daughterHelixClosestToMotherPosLast.getD0() * daughterHelixClosestToMotherPosLast.getD0() +
           daughterHelixClosestToMotherPosLast.getZ0() * daughterHelixClosestToMotherPosLast.getZ0()) <
          m_precutDistance * m_precutDistance);
}

refitRecoTrackAfterReassign()

bool refitRecoTrackAfterReassign ( RecoTrack *  recoTrackMotherRefit, RecoTrack *  recoTrackDaughterRefit, const RecoTrack *  recoTrackMother, const RecoTrack *  recoTrackDaughter  )

private

Try to fit new RecoTracks after hit reassignment.

Parameters

recoTrackMotherRefit	mother track after reassignment
recoTrackDaughterRefit	daughter track after reassignment
recoTrackMother	initial mother track
recoTrackDaughter	initial daughter track

Returns

true if the fit was successful and the result for two tracks improved

Definition at line 311 of file KinkFitter.cc.

{
 
  // initialize fitter
  TrackFitter trackFitterDAF;
  // fit the new tracks
  trackFitterDAF.fit(*recoTrackMotherRefit);
  trackFitterDAF.fit(*recoTrackDaughterRefit);
 
  if (!recoTrackMotherRefit->wasFitSuccessful() || !recoTrackDaughterRefit->wasFitSuccessful()) {
    B2DEBUG(29, "Refit of the tracks with reassigned hits failed ");
    return false;
  }
 
  const genfit::FitStatus* motherTrackFitStatus = recoTrackMother->getTrackFitStatus();
  const genfit::FitStatus* daughterTrackFitStatus = recoTrackDaughter->getTrackFitStatus();
  const double chi2MotherInit = motherTrackFitStatus->getChi2();
  const double chi2DaughterInit = daughterTrackFitStatus->getChi2();
  const double ndfMotherInit = motherTrackFitStatus->getNdf();
  const double ndfDaughterInit = daughterTrackFitStatus->getNdf();
  B2DEBUG(29, "Initial mother fit result, p-value: " << motherTrackFitStatus->getPVal() << ", chi2: " <<
          chi2MotherInit << ", ndf: " << ndfMotherInit);
  B2DEBUG(29, "Initial daughter fit result, p-value: " << daughterTrackFitStatus->getPVal() << ", chi2: " <<
          chi2DaughterInit << ", ndf: " << ndfDaughterInit);
 
  const genfit::FitStatus* motherNewTrackFitStatus = recoTrackMotherRefit->getTrackFitStatus();
  const genfit::FitStatus* daughterNewTrackFitStatus = recoTrackDaughterRefit->getTrackFitStatus();
  const double chi2Mother = motherNewTrackFitStatus->getChi2();
  const double chi2Daughter = daughterNewTrackFitStatus->getChi2();
  const double ndfMother = motherNewTrackFitStatus->getNdf();
  const double ndfDaughter = daughterNewTrackFitStatus->getNdf();
  B2DEBUG(29, "New mother fit result, p-value: " << motherNewTrackFitStatus->getPVal() << ", chi2: " <<
          chi2Mother << ", ndf: " << ndfMother);
  B2DEBUG(29, "New daughter fit result, p-value: " << daughterNewTrackFitStatus->getPVal() << ", chi2: " <<
          chi2Daughter << ", ndf: " << ndfDaughter);
 
  // Fit is assumed to be successful if it improves sum of the chi2 divided by sum of ndf of two tracks
  if ((chi2Mother + chi2Daughter) / (ndfMother + ndfDaughter) < (chi2MotherInit + chi2DaughterInit) /
      (ndfMotherInit + ndfDaughterInit))
    return true;
  else
    return false;
 
}

setFitterMode()

void setFitterMode

(

const unsigned char

fitterMode

)

set kink fitter mode.

Parameters

fitterMode

from 0 to 15 in binary: 1st bit: reassign hits (1 is On, 0 is Off) 2nd bit: flip tracks with close end points (1 is On, 0 is Off) 3rd bit: fit both tracks as one (1 is On, 0 is Off) 4th bit: track splitting (1 is On, 0 is Off)

Definition at line 68 of file KinkFitter.cc.

{
  if (not(fitterMode <= 15)) {
    B2FATAL("Invalid fitter mode!");
  } else {
    m_kinkFitterMode = fitterMode;
    // filling bits of the fitter mode
    m_kinkFitterModeHitsReassignment = fitterMode & 0b0001;
    m_kinkFitterModeFlipAndRefit = fitterMode & 0b0010;
    m_kinkFitterModeCombineAndFit = fitterMode & 0b0100;
    m_kinkFitterModeSplitTrack = fitterMode & 0b1000;
  }
}

splitRecoTrack()

bool splitRecoTrack ( const RecoTrack *  recoTrackSplit, short &  recoTrackIndexMother, short &  recoTrackIndexDaughter  )

private

Split track into two based on |chi2/ndf - 1|.

For the best split position search, use binary search.

Parameters

recoTrackSplit	RecoTrack to be split
recoTrackIndexMother	index of the created mother RecoTrack in m_copiedRecoTracks
recoTrackIndexDaughter	index of the created daughter RecoTrack in m_copiedRecoTracks

Returns

true if the splitting is successful, false otherwise

Definition at line 612 of file KinkFitter.cc.

{
  // initialize fitter
  TrackFitter trackFitterDAF;
 
  // initial default value for the |chi2/ndf - 1|
  constexpr double defaultChi2NdfRatio = std::numeric_limits<double>::max();
 
  // enum for the splitting points used in the best splitting point search
  enum {
    c_outerEdge, c_middlePoint, c_innerEdge
  };
 
  B2DEBUG(29, "Start splitting");
 
  // fit result of the initial track to be split
  const genfit::FitStatus* splitTrackFitStatus = recoTrackSplit->getTrackFitStatus();
  const double ndfSplit = splitTrackFitStatus->getNdf();
  const double chi2Split = splitTrackFitStatus->getChi2();
  const double chi2NdfRatioSplit = fabs(chi2Split / ndfSplit - 1);
  B2DEBUG(29, "Initial fit result, p-value: " <<  splitTrackFitStatus->getPVal() << ", chi2: " <<
          chi2Split << ", ndf: " << ndfSplit);
 
  // number of SVD+CDC hits of the initial track to be split
  const unsigned int numberSVDCDCHitsSplit = recoTrackSplit->getNumberOfCDCHits() + recoTrackSplit->getNumberOfSVDHits();
 
  // Three initial points for the binary search. Since the main target are tracks with mixed fractions of hits from two
  // wrongly combined tracks between 30-70%, the points are chosen between 20-80%.
  // The first point is for 80% (20%) of hits assigned to mother (daughter), the second one is for middle 50% (50%),
  // and the last one is for 20% (80%), respectively.
  // Due to the logic of the copyRecoTrackAndSplit function, the hits position in splitHitNumber
  // is counted from the end of the track to split (the index 0 is for c_outerEdge, index 1 is for c_middlePoint,
  // and index 2 is for c_innerEdge).
  std::array<unsigned int, 3> splitHitNumber = {
    static_cast<unsigned int>(numberSVDCDCHitsSplit * 0.2),
    static_cast<unsigned int>(numberSVDCDCHitsSplit * 0.5),
    static_cast<unsigned int>(numberSVDCDCHitsSplit * 0.8)
  };
 
  // initialization of arrays of the created mother and daughter RecoTracks and |chi2/ndf - 1| for them
  // chi2/ndf ratios for resulting tracks are calculated as (chi2_mother + chi2_daughter)/(ndf_mother + ndf_daughter)
  std::array<RecoTrack*, 3> recoTrackMother = {nullptr, nullptr, nullptr};
  std::array<RecoTrack*, 3> recoTrackDaughter = {nullptr, nullptr, nullptr};
  std::array<double, 3> chi2NdfRatio = {defaultChi2NdfRatio, defaultChi2NdfRatio, defaultChi2NdfRatio};
 
  // initial splitting and fitting of the resulting tracks
  for (size_t splitHitNumberIndex = 0; splitHitNumberIndex < splitHitNumber.size(); ++splitHitNumberIndex) {
    // if there is not enough hits for fit of at least one resulting track, continue
    if (splitHitNumber[splitHitNumberIndex] < 5 ||
        (numberSVDCDCHitsSplit - splitHitNumber[splitHitNumberIndex]) < 5) continue;
 
    // Splitting is based only on the CDC hits, so all VXD hits are assigned to mother. If the splitting number exceeds
    // number of CDC hits, use the latter one.
    if (splitHitNumber[splitHitNumberIndex] > recoTrackSplit->getNumberOfCDCHits())
      splitHitNumber[splitHitNumberIndex] = recoTrackSplit->getNumberOfCDCHits();
 
    B2DEBUG(29, "Initial splitting hit number for point " << splitHitNumberIndex << ": " <<
            splitHitNumber[splitHitNumberIndex]);
 
    // creation of split RecoTracks
    recoTrackMother[splitHitNumberIndex] = copyRecoTrackAndSplit(recoTrackSplit, true, splitHitNumber[splitHitNumberIndex]);
    recoTrackDaughter[splitHitNumberIndex] = copyRecoTrackAndSplit(recoTrackSplit, false, splitHitNumber[splitHitNumberIndex]);
 
    // fit of the RecoTracks
    trackFitterDAF.fit(*(recoTrackMother[splitHitNumberIndex]));
    trackFitterDAF.fit(*(recoTrackDaughter[splitHitNumberIndex]));
 
    // if the fit failed, skip this track pair
    if (!recoTrackMother[splitHitNumberIndex]->wasFitSuccessful() ||
        !recoTrackDaughter[splitHitNumberIndex]->wasFitSuccessful()) continue;
 
    // fit result of the split track pair
    const genfit::FitStatus* motherNewTrackFitStatus = recoTrackMother[splitHitNumberIndex]->getTrackFitStatus();
    const genfit::FitStatus* daughterNewTrackFitStatus = recoTrackDaughter[splitHitNumberIndex]->getTrackFitStatus();
    const double ndfMotherTmp = motherNewTrackFitStatus->getNdf();
    const double ndfDaughterTmp = daughterNewTrackFitStatus->getNdf();
    const double chi2MotherTmp = motherNewTrackFitStatus->getChi2();
    const double chi2DaughterTmp = daughterNewTrackFitStatus->getChi2();
    const double chi2NdfRatioTmp = fabs((chi2MotherTmp + chi2DaughterTmp) / (ndfMotherTmp + ndfDaughterTmp) - 1);
    B2DEBUG(29, "Mother fit result for initial point " << splitHitNumberIndex << ", p-value: " <<
            motherNewTrackFitStatus->getPVal() << ", chi2: " << chi2MotherTmp << ", ndf: " << ndfMotherTmp);
    B2DEBUG(29, "Daughter fit result for initial point " << splitHitNumberIndex << ", p-value: " <<
            daughterNewTrackFitStatus->getPVal() << ", chi2: " << chi2DaughterTmp << ", ndf: " << ndfDaughterTmp);
    B2DEBUG(29, "|Chi2/NDF - 1| for initial point " << splitHitNumberIndex << ": " <<
            chi2NdfRatioTmp << ", and for initial track: " << chi2NdfRatioSplit);
 
    // if the resulting |chi2/ndf - 1| is bigger than one for initial track, skip this track pair
    if (chi2NdfRatioTmp > chi2NdfRatioSplit) continue;
 
    chi2NdfRatio[splitHitNumberIndex] = chi2NdfRatioTmp;
  }
 
  // Find the point with minimal difference of chi2/ndf ratio with unity using binary search.
  // Compare the |chi2/ndf - 1| result for edge points and choose the interval between middle and the better edge.
  // Find a middle point in the chosen interval, split the track, and fit the result.
  // Return to the comparison.
  // In a special case when both edges fit failed (the values of |chi2/ndf - 1| are equal),
  // check if the middle point is fitted. If it is, use it as a result. If it is not, return false.
  // Try five iterations (usually, enough to converge). If the middle point equals one of the edge, the result has
  // converged earlier. Break the loop and proceed the execution of the function.
  for (unsigned int searchIteration = 0; searchIteration < 5; ++searchIteration) {
    B2DEBUG(29, "Splitting iteration " << searchIteration << "; hit numbers: " << splitHitNumber[c_outerEdge] << ", " <<
            splitHitNumber[c_middlePoint] << ", " << splitHitNumber[c_innerEdge]);
 
    if (chi2NdfRatio[c_outerEdge] < chi2NdfRatio[c_innerEdge]) { // choose outer interval
 
      // replace the inner edge with the middle point
      splitHitNumber[c_innerEdge] = splitHitNumber[c_middlePoint];
      chi2NdfRatio[c_innerEdge] = chi2NdfRatio[c_middlePoint];
      recoTrackMother[c_innerEdge] = recoTrackMother[c_middlePoint];
      recoTrackDaughter[c_innerEdge] = recoTrackDaughter[c_middlePoint];
      // find a new middle point in the outer interval, and fill the other arrays with the default values
      splitHitNumber[c_middlePoint] = (splitHitNumber[c_innerEdge] + splitHitNumber[c_outerEdge]) / 2;
      chi2NdfRatio[c_middlePoint] = defaultChi2NdfRatio;
      recoTrackMother[c_middlePoint] = nullptr;
      recoTrackDaughter[c_middlePoint] = nullptr;
 
      B2DEBUG(29, "Outer interval is chosen; hit numbers: " << splitHitNumber[c_outerEdge] << ", " <<
              splitHitNumber[c_middlePoint] << ", " << splitHitNumber[c_innerEdge]);
 
    } else if (chi2NdfRatio[c_outerEdge] > chi2NdfRatio[c_innerEdge]) { // choose inner interval
 
      // replace the outer edge with the middle point
      splitHitNumber[c_outerEdge] = splitHitNumber[c_middlePoint];
      chi2NdfRatio[c_outerEdge] = chi2NdfRatio[c_middlePoint];
      recoTrackMother[c_outerEdge] = recoTrackMother[c_middlePoint];
      recoTrackDaughter[c_outerEdge] = recoTrackDaughter[c_middlePoint];
      // find a new middle point in the inner interval, and fill the other arrays with the default values
      splitHitNumber[c_middlePoint] = (splitHitNumber[c_innerEdge] + splitHitNumber[c_outerEdge]) / 2;
      chi2NdfRatio[c_middlePoint] = defaultChi2NdfRatio;
      recoTrackMother[c_middlePoint] = nullptr;
      recoTrackDaughter[c_middlePoint] = nullptr;
 
      B2DEBUG(29, "Inner interval is chosen; hit numbers: " << splitHitNumber[c_outerEdge] << ", " <<
              splitHitNumber[c_middlePoint] << ", " << splitHitNumber[c_innerEdge]);
 
    } else if (chi2NdfRatio[c_middlePoint] == defaultChi2NdfRatio) { // if none of the point was fitted successfully, return false
 
      B2DEBUG(29, "Track splitting failed");
      return false;
 
    } else { // if both edges fits have the same results (their fit failed), return the middle one
 
      B2DEBUG(29, "Middle point is chosen; hit numbers: " << splitHitNumber[c_outerEdge] << ", " <<
              splitHitNumber[c_middlePoint] << ", " << splitHitNumber[c_innerEdge]);
      // return through references indexes of the mother and daughter tracks
      recoTrackIndexMother = recoTrackMother[c_middlePoint]->getArrayIndex();
      recoTrackIndexDaughter = recoTrackDaughter[c_middlePoint]->getArrayIndex();
      return true;
 
    }
 
    // stop the iterations if the all possibilities were tried
    if (splitHitNumber[c_middlePoint] == splitHitNumber[c_outerEdge] ||
        splitHitNumber[c_middlePoint] == splitHitNumber[c_innerEdge]) break;
 
    // if there is not enough hits for fit of at least one resulting track, continue
    if (splitHitNumber[c_middlePoint] < 5 ||
        (numberSVDCDCHitsSplit - splitHitNumber[c_middlePoint]) < 5) continue;
 
    // creation of split RecoTracks in the new middle point
    recoTrackMother[c_middlePoint] = copyRecoTrackAndSplit(recoTrackSplit, true, splitHitNumber[c_middlePoint]);
    recoTrackDaughter[c_middlePoint] = copyRecoTrackAndSplit(recoTrackSplit, false, splitHitNumber[c_middlePoint]);
 
    // fit of the mother and daughter RecoTracks in the new middle point
    trackFitterDAF.fit(*(recoTrackMother[c_middlePoint]));
    trackFitterDAF.fit(*(recoTrackDaughter[c_middlePoint]));
 
    // if the fit failed, skip this track pair
    if (!recoTrackMother[c_middlePoint]->wasFitSuccessful() ||
        !recoTrackDaughter[c_middlePoint]->wasFitSuccessful()) continue;
 
    // fit result of the split track pair
    const genfit::FitStatus* motherNewTrackFitStatus = recoTrackMother[c_middlePoint]->getTrackFitStatus();
    const genfit::FitStatus* daughterNewTrackFitStatus = recoTrackDaughter[c_middlePoint]->getTrackFitStatus();
    const double ndfMotherTmp = motherNewTrackFitStatus->getNdf();
    const double ndfDaughterTmp = daughterNewTrackFitStatus->getNdf();
    const double chi2MotherTmp = motherNewTrackFitStatus->getChi2();
    const double chi2DaughterTmp = daughterNewTrackFitStatus->getChi2();
    const double chi2NdfRatioTmp = fabs((chi2MotherTmp + chi2DaughterTmp) / (ndfMotherTmp + ndfDaughterTmp) - 1);
 
    B2DEBUG(29, "Mother fit result for approach " << searchIteration << ", p-value: " <<
            motherNewTrackFitStatus->getPVal() << ", chi2: " << chi2MotherTmp << ", ndf: " << ndfMotherTmp);
    B2DEBUG(29, "Daughter fit result for approach " << searchIteration << ", p-value: " <<
            daughterNewTrackFitStatus->getPVal() << ", chi2: " << chi2DaughterTmp << ", ndf: " << ndfDaughterTmp);
    B2DEBUG(29, "Chi2/NDF for approach " << searchIteration << ": " <<
            chi2NdfRatioTmp << ", and for initial track: " << chi2NdfRatioSplit);
 
    // if the resulting chi2/ndf ratio is bigger than one for initial track, skip this track pair
    if (chi2NdfRatioTmp > chi2NdfRatioSplit) continue;
 
    chi2NdfRatio[c_middlePoint] = chi2NdfRatioTmp;
  }
 
  // find the minimal difference of chi2/ndf ratio with unity among left points
  std::array<double, 3>::iterator minChi2NdfRatioResult = std::min_element(chi2NdfRatio.begin(), chi2NdfRatio.end());
  const size_t minIndex = std::distance(chi2NdfRatio.begin(), minChi2NdfRatioResult);
 
 
  B2DEBUG(29, "End of splitting with |chi2/ndf - 1| for point 1:" << chi2NdfRatio[c_outerEdge] <<
          "; 2:" << chi2NdfRatio[c_middlePoint] << "; 3:" << chi2NdfRatio[c_innerEdge]);
  B2DEBUG(29, "Hit numbers for point 1:" << splitHitNumber[c_outerEdge] <<
          "; 2:" << splitHitNumber[c_middlePoint] << "; 3:" << splitHitNumber[c_innerEdge]);
  B2DEBUG(29, "Min index: " << minIndex);
 
  // return through references indexes of the mother and daughter tracks for the best splitting point
  recoTrackIndexMother = recoTrackMother[minIndex]->getArrayIndex();
  recoTrackIndexDaughter = recoTrackDaughter[minIndex]->getArrayIndex();
 
  return true;
 
}

vertexFitWithRecoTracks()

bool vertexFitWithRecoTracks ( RecoTrack *  recoTrackMother, RecoTrack *  recoTrackDaughter, unsigned int &  reassignHitStatus, ROOT::Math::XYZVector &  vertexPos, double &  distance, ROOT::Math::XYZVector  vertexPosSeed = ROOT::Math::XYZVector(0, 0, 0)  )

private

Fit kink vertex using RecoTrack's as inputs.

Fit kink vertex using RecoTracks as inputs.

Return true (false) if the vertex fit has done well (failed). If the vertex is inside one of the RecoTracks, bits in reassignHitStatus are set.

Parameters

recoTrackMother	RecoTrack of mother
recoTrackDaughter	RecoTrack of daughter
reassignHitStatus	store a result of this function. if the daughter(mother) track has hits inside the Kink vertex position, 0x1(0x2) bit is set.
vertexPos	store a result of this function. The fitted vertex position is stored.
distance	store a distance between tracks at the decay vertex
vertexPosSeed	a seed of the vertex position

Returns

true if the vertex fit succeed; false if one of the tracks is not fitted (should not happen), kFit finishes with error, chi^2 of the fit is too large, or extrapolation of one of the tracks to decay vertex fails

Calculates distance at the fitted vertex. Checks if the reassignment of the hits is required.

Definition at line 1280 of file KinkFitter.cc.

{
 
  // make a clone, not use the reference so that the genfit::Track and its TrackReps will not be altered.
  genfit::AbsTrackRep* motherRepresentation = recoTrackMother->getCardinalRepresentation();
  if ((motherRepresentation == nullptr) || !(recoTrackMother->wasFitSuccessful(motherRepresentation))) {
    B2ERROR("Cardinal representation is not available for track. Should never happen, but I can continue safely anyway.");
    return false;
  }
 
  const double motherMass = TDatabasePDG::Instance()->GetParticle(motherRepresentation->getPDG())->Mass();
  const double motherCharge = recoTrackMother->getTrackFitStatus()->getCharge();
 
  // make a clone, not use the reference so that the genfit::Track and its TrackReps will not be altered.
  genfit::AbsTrackRep* daughterRepresentation = recoTrackDaughter->getCardinalRepresentation();
  if ((daughterRepresentation == nullptr) || !(recoTrackDaughter->wasFitSuccessful(daughterRepresentation))) {
    B2ERROR("Cardinal representation is not available for track. Should never happen, but I can continue safely anyway.");
    return false;
  }
 
  const double daughterMass = TDatabasePDG::Instance()->GetParticle(daughterRepresentation->getPDG())->Mass();
  const double daughterCharge = recoTrackDaughter->getTrackFitStatus()->getCharge();
 
  // make a clone, not use the reference so that the genfit::MeasuredStateOnPlane and its TrackReps will not be altered.
  genfit::MeasuredStateOnPlane stMother = recoTrackMother->getMeasuredStateOnPlaneFromLastHit(motherRepresentation);
  genfit::MeasuredStateOnPlane stDaughter = recoTrackDaughter->getMeasuredStateOnPlaneClosestTo(vertexPosSeed,
                                            daughterRepresentation);
  m_stMotherBuffer = stMother;
  m_stDaughterBuffer = stDaughter;
 
  // initialize KVertexFitter
  analysis::VertexFitKFit kvf;
 
  // set magnetic field at the seed position
  const double Bz = BFieldManager::getFieldInTesla(vertexPosSeed).Z();
  kvf.setMagneticField(Bz);
 
  // set seed position
  HepPoint3D vertexPosSeedHepPoint(vertexPosSeed.X(), vertexPosSeed.Y(), vertexPosSeed.Z());
  kvf.setInitialVertex(vertexPosSeedHepPoint);
 
  // add mother state for the fit
  TVector3 motherPosition, motherMomentum;
  TMatrixDSym motherCovMatrix6(6, 6);
  stMother.getPosMomCov(motherPosition, motherMomentum, motherCovMatrix6);
 
  const double motherEnergy = sqrt(motherMomentum.Mag2() + motherMass * motherMass);
  ROOT::Math::PxPyPzEVector motherFourMomentum(motherMomentum.X(), motherMomentum.Y(),
                                               motherMomentum.Z(), motherEnergy);
 
  TMatrixDSym motherErrMatrixForKFit(7);
  errMatrixForKFit(motherFourMomentum, motherCovMatrix6, motherErrMatrixForKFit);
 
  kvf.addTrack(ROOTToCLHEP::getHepLorentzVector(motherFourMomentum),
               ROOTToCLHEP::getPoint3D(ROOT::Math::XYZVector(motherPosition)),
               ROOTToCLHEP::getHepSymMatrix(motherErrMatrixForKFit),
               motherCharge);
 
  // add daughter state for the fit
  TVector3 daughterPosition, daughterMomentum;
  TMatrixDSym daughterCovMatrix6(6, 6);
  stDaughter.getPosMomCov(daughterPosition, daughterMomentum, daughterCovMatrix6);
 
  const double daughterEnergy = sqrt(daughterMomentum.Mag2() + daughterMass * daughterMass);
  ROOT::Math::PxPyPzEVector daughterFourMomentum(daughterMomentum.X(), daughterMomentum.Y(),
                                                 daughterMomentum.Z(), daughterEnergy);
 
  TMatrixDSym daughterErrMatrixForKFit(7);
  errMatrixForKFit(daughterFourMomentum, daughterCovMatrix6, daughterErrMatrixForKFit);
 
  kvf.addTrack(ROOTToCLHEP::getHepLorentzVector(daughterFourMomentum),
               ROOTToCLHEP::getPoint3D(ROOT::Math::XYZVector(daughterPosition)),
               ROOTToCLHEP::getHepSymMatrix(daughterErrMatrixForKFit),
               daughterCharge);
 
  // do the fit
  int err = kvf.doFit();
  if (err != 0) {
    B2DEBUG(29, "Vertex fit finished with error");
    return false;
  }
 
  // test the chi2 cut
  if (kvf.getCHIsq() > m_vertexChi2Cut) {
    B2DEBUG(29, "Chi^2  of vertex fit is too large: " << kvf.getCHIsq());
    return false;
  }
 
  // get the fitted vertex
  HepPoint3D vertexPosHepPoint = kvf.getVertex();
  vertexPos = ROOT::Math::XYZVector(vertexPosHepPoint.x(), vertexPosHepPoint.y(), vertexPosHepPoint.z());
 
  // extrapolate the mother and the daughter states to the fitted vertex and get the status of the hit reassignment
  if (!extrapolateToVertex(stMother, stDaughter, vertexPos, reassignHitStatus)) {
    B2DEBUG(29, "Failed to extrapolate one of the tracks to the fitted kink vertex");
    return false;
  }
 
  // prepare mother Helix at the fitted vertex to calculate the distance between tracks
  const ROOT::Math::XYZVector motherPos = ROOT::Math::XYZVector(stMother.getPos());
  const ROOT::Math::XYZVector motherMom = ROOT::Math::XYZVector(stMother.getMom());
  Helix motherHelix(motherPos,
                    motherMom,
                    static_cast<short>(recoTrackMother->getTrackFitStatus()->getCharge()),
                    Bz);
  motherHelix.passiveMoveBy(vertexPos);
 
  // prepare daughter Helix at the fitted vertex to calculate the distance between tracks
  const ROOT::Math::XYZVector daughterPos = ROOT::Math::XYZVector(stDaughter.getPos());
  const ROOT::Math::XYZVector daughterMom = ROOT::Math::XYZVector(stDaughter.getMom());
  Helix daughterHelix(daughterPos,
                      daughterMom,
                      static_cast<short>(recoTrackDaughter->getTrackFitStatus()->getCharge()),
                      Bz);
  daughterHelix.passiveMoveBy(vertexPos);
 
  // calculate the distance between tracks as squared sum of their impact parameters at the fitted vertex
  distance = sqrt(daughterHelix.getD0() * daughterHelix.getD0() + daughterHelix.getZ0() * daughterHelix.getZ0() +
                  motherHelix.getD0() * motherHelix.getD0() + motherHelix.getZ0() * motherHelix.getZ0());
 
  // save mother and daughter states to the buffer
  m_stMotherBuffer = stMother;
  m_stDaughterBuffer = stDaughter;
 
  B2DEBUG(29, "Vertex fit with chi2 " << kvf.getCHIsq() << " and distance between tracks " << distance);
  return true;
}

Member Data Documentation

m_copiedRecoTracks

StoreArray<RecoTrack> m_copiedRecoTracks

private

RecoTrack used to refit tracks.

Definition at line 282 of file KinkFitter.h.

m_daughterKinkRecoTrackCache

RecoTrack* m_daughterKinkRecoTrackCache

private

cache for the RecoTrack of daughter used to find the best vertex

Definition at line 286 of file KinkFitter.h.

m_kinkFitterMode

unsigned char m_kinkFitterMode

private

fitter mode from 0 to 15 written in bits:

1st bit: reassign hits (1 is On, 0 is Off) 2nd bit: flip tracks with close end points (1 is On, 0 is Off) 3rd bit: fit both tracks as one (1 is On, 0 is Off) 4th bit: combined track candidate splitting (1 is On, 0 is Off)

Definition at line 271 of file KinkFitter.h.

m_kinkFitterModeCombineAndFit

bool m_kinkFitterModeCombineAndFit

private

fitter mode 3rd bit

Definition at line 278 of file KinkFitter.h.

m_kinkFitterModeFlipAndRefit

bool m_kinkFitterModeFlipAndRefit

private

fitter mode 2nd bit

Definition at line 277 of file KinkFitter.h.

m_kinkFitterModeHitsReassignment

bool m_kinkFitterModeHitsReassignment

private

fitter mode 1st bit

Definition at line 276 of file KinkFitter.h.

m_kinkFitterModeSplitTrack

bool m_kinkFitterModeSplitTrack

private

fitter mode 4th bit

Definition at line 279 of file KinkFitter.h.

m_kinks

StoreArray<Kink> m_kinks

private

Kink (output).

Definition at line 262 of file KinkFitter.h.

m_motherKinkRecoTrackCache

RecoTrack* m_motherKinkRecoTrackCache

private

cache for the RecoTrack of mother used to find the best vertex

Definition at line 285 of file KinkFitter.h.

m_precutDistance

double m_precutDistance

private

Preselection cut on distance between ending points of two tracks used in prefilter.

here it is needed in isRefitImproveFilter6 function

Definition at line 267 of file KinkFitter.h.

m_recoTracks

StoreArray<RecoTrack> m_recoTracks

private

RecoTrack (input)

Definition at line 258 of file KinkFitter.h.

m_recoTracksName

std::string m_recoTracksName

private

RecoTrackColName (input).

Definition at line 257 of file KinkFitter.h.

m_stDaughterBuffer

genfit::MeasuredStateOnPlane m_stDaughterBuffer

private

buffer for the MeasuredStateOnPlane of daughter obtained in the vertex fit

Definition at line 284 of file KinkFitter.h.

m_stMotherBuffer

genfit::MeasuredStateOnPlane m_stMotherBuffer

private

buffer for the MeasuredStateOnPlane of mother obtained in the vertex fit

Definition at line 283 of file KinkFitter.h.

m_trackFitResults

StoreArray<TrackFitResult> m_trackFitResults

private

TrackFitResult (output).

Definition at line 261 of file KinkFitter.h.

m_vertexChi2Cut

double m_vertexChi2Cut

private

Chi2 cut.

Definition at line 266 of file KinkFitter.h.

m_vertexDistanceCut

double m_vertexDistanceCut

private

cut on the distance at the found vertex.

Definition at line 265 of file KinkFitter.h.

---

The documentation for this class was generated from the following files:

• tracking/kinkFinding/fitter/include/KinkFitter.h
• tracking/kinkFinding/fitter/src/KinkFitter.cc