V0Fitter.cc

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
#include <tracking/v0Finding/fitter/V0Fitter.h>
 
#include <framework/logging/Logger.h>
#include <framework/datastore/StoreArray.h>
#include <framework/geometry/VectorUtil.h>
#include <framework/geometry/BFieldManager.h>
#include <tracking/v0Finding/dataobjects/VertexVector.h>
#include <tracking/dataobjects/RecoTrack.h>
#include <tracking/trackFitting/fitter/base/TrackFitter.h>
#include <tracking/trackFitting/trackBuilder/factories/TrackBuilder.h>
 
#include <mdst/dataobjects/HitPatternVXD.h>
#include <mdst/dataobjects/HitPatternCDC.h>
#include <mdst/dataobjects/Track.h>
#include <mdst/dataobjects/TrackFitResult.h>
 
#include <genfit/Track.h>
#include <genfit/TrackPoint.h>
#include <genfit/MeasuredStateOnPlane.h>
#include <genfit/GFRaveVertexFactory.h>
#include <genfit/GFRaveVertex.h>
#include <genfit/FieldManager.h>
#include <genfit/MaterialEffects.h>
#include <genfit/FitStatus.h>
#include <genfit/KalmanFitStatus.h>
#include <genfit/KalmanFitterInfo.h>
 
#include <framework/utilities/IOIntercept.h>
 
using namespace Belle2;
 
V0Fitter::V0Fitter(const std::string& trackFitResultsName, const std::string& v0sName,
                   const std::string& v0ValidationVerticesName, const std::string& recoTracksName,
                   const std::string& copiedRecoTracksName, bool enableValidation)
  : m_validation(enableValidation), m_recoTracksName(recoTracksName), m_v0FitterMode(1), m_forcestore(false),
    m_useOnlyOneSVDHitPair(true)
{
  m_trackFitResults.isRequired(trackFitResultsName);
  m_v0s.registerInDataStore(v0sName, DataStore::c_WriteOut | DataStore::c_ErrorIfAlreadyRegistered);
  //Relation to RecoTracks from Tracks is already tested at the module level.
 
  if (m_validation) {
    B2DEBUG(24, "Register DataStore for validation.");
    m_validationV0s.registerInDataStore(v0ValidationVerticesName, DataStore::c_ErrorIfAlreadyRegistered);
    m_v0s.registerRelationTo(m_validationV0s);
  }

  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());
}
 
void V0Fitter::setFitterMode(int fitterMode)
{
  if (not(0 <= fitterMode && fitterMode <= 2)) {
    B2FATAL("Invalid fitter mode!");
  } else {
    m_v0FitterMode = fitterMode;
    if (fitterMode == 0)
      m_forcestore = true;
    else
      m_forcestore = false;
    if (fitterMode == 2)
      m_useOnlyOneSVDHitPair = false;
    else
      m_useOnlyOneSVDHitPair = true;
  }
}
 
void V0Fitter::initializeCuts(double beamPipeRadius,
                              double vertexChi2CutOutside,
                              std::tuple<double, double> invMassRangeKshort,
                              std::tuple<double, double> invMassRangeLambda,
                              std::tuple<double, double> invMassRangePhoton)
{
  m_beamPipeRadius = beamPipeRadius;
  m_vertexChi2CutOutside = vertexChi2CutOutside;
  m_invMassRangeKshort = invMassRangeKshort;
  m_invMassRangeLambda = invMassRangeLambda;
  m_invMassRangePhoton = invMassRangePhoton;
}
 
 
// cppcheck-suppress[constParameterReference] ; genfit::GFRaveVertexFactory::findVertices takes non-const Track pointers
bool V0Fitter::fitGFRaveVertex(genfit::Track& trackPlus, genfit::Track& trackMinus, genfit::GFRaveVertex& vertex)
{
  VertexVector vertexVector;
  std::vector<genfit::Track*> trackPair {&trackPlus, &trackMinus};
 
  try {
    IOIntercept::OutputToLogMessages
    logCapture("V0Fitter GFRaveVertexFactory", LogConfig::c_Debug, LogConfig::c_Debug);
    logCapture.start();
 
    genfit::GFRaveVertexFactory vertexFactory;
    vertexFactory.findVertices(&vertexVector.v, trackPair);
  } catch (...) {
    B2ERROR("Exception during vertex fit.");
    return false;
  }
 
  if (vertexVector.size() != 1) {
    B2DEBUG(21, "Vertex fit failed. Size of vertexVector not 1, but: " << vertexVector.size());
    return false;
  }
 
  if ((*vertexVector[0]).getNTracks() != 2) {
    B2DEBUG(20, "Wrong number of tracks in vertex.");
    return false;
  }
 
  vertex = *vertexVector[0];
  return true;
}

bool V0Fitter::extrapolateToVertex(genfit::MeasuredStateOnPlane& stPlus, genfit::MeasuredStateOnPlane& stMinus,
                                   const ROOT::Math::XYZVector& vertexPosition, unsigned int& hasInnerHitStatus)
{
  hasInnerHitStatus = 0;
 
  try {
    const double extralengthPlus  =  stPlus.extrapolateToPoint(XYZToTVector(vertexPosition));
    const double extralengthMinus = stMinus.extrapolateToPoint(XYZToTVector(vertexPosition));
    if (extralengthPlus  > 0) hasInnerHitStatus |= 0x1; 
    if (extralengthMinus > 0) hasInnerHitStatus |= 0x2; 
    B2DEBUG(22, "extralengthPlus=" << extralengthPlus << ", extralengthMinus=" << extralengthMinus);
  } catch (...) {
    B2DEBUG(22, "Could not extrapolate track to vertex.");
    return false;
  }
  return true;
}
 
TrackFitResult* V0Fitter::buildTrackFitResult(const genfit::Track& track, const RecoTrack* recoTrack,
                                              const genfit::MeasuredStateOnPlane& msop, const double Bz,
                                              const Const::ParticleType& trackHypothesis,
                                              const int sharedInnermostCluster)
{
  const uint64_t hitPatternCDCInitializer = TrackBuilder::getHitPatternCDCInitializer(*recoTrack);
  uint32_t hitPatternVXDInitializer = TrackBuilder::getHitPatternVXDInitializer(*recoTrack);
 
  // If the innermost hit is shared among V0 daughters, assign flag in the infoLayer.
  if (sharedInnermostCluster > 0 && sharedInnermostCluster < 4) {
    HitPatternVXD hitPatternVXD_forflag = HitPatternVXD(hitPatternVXDInitializer);
    hitPatternVXD_forflag.setInnermostHitShareStatus(sharedInnermostCluster);
    hitPatternVXDInitializer = hitPatternVXD_forflag.getInteger();
  }
 
  TrackFitResult* v0TrackFitResult
    = m_trackFitResults.appendNew(ROOT::Math::XYZVector(msop.getPos()), ROOT::Math::XYZVector(msop.getMom()),
                                  msop.get6DCov(), msop.getCharge(),
                                  trackHypothesis,
                                  track.getFitStatus()->getPVal(),
                                  Bz, hitPatternCDCInitializer, hitPatternVXDInitializer, track.getFitStatus()->getNdf());
  return v0TrackFitResult;
}
 
std::pair<Const::ParticleType, Const::ParticleType> V0Fitter::getTrackHypotheses(const Const::ParticleType& v0Hypothesis)
{
  if (v0Hypothesis == Const::Kshort) {
    return std::make_pair(Const::pion, Const::pion);
  } else if (v0Hypothesis == Const::photon) {
    return std::make_pair(Const::electron, Const::electron);
  } else if (v0Hypothesis == Const::Lambda) {
    return std::make_pair(Const::proton, Const::pion);
  } else if (v0Hypothesis == Const::antiLambda) {
    return std::make_pair(Const::pion, Const::proton);
  }
  B2FATAL("Given V0Hypothesis not available.");
  return std::make_pair(Const::invalidParticle, Const::invalidParticle); // return something to avoid triggering cppcheck
}

bool V0Fitter::fitAndStore(const Track* trackPlus, const Track* trackMinus,
                           const Const::ParticleType& v0Hypothesis,
                           bool& isForceStored, bool& isHitRemoved)
{
  isForceStored = false;
  isHitRemoved = false;

  RecoTrack* recoTrackPlus  =  trackPlus->getRelated<RecoTrack>(m_recoTracksName);
  RecoTrack* recoTrackMinus = trackMinus->getRelated<RecoTrack>(m_recoTracksName);

  unsigned int hasInnerHitStatus = 0;

  ROOT::Math::XYZVector vertexPos(0, 0, 0);

  if (not vertexFitWithRecoTracks(trackPlus, trackMinus, recoTrackPlus, recoTrackMinus, v0Hypothesis,
                                  hasInnerHitStatus, vertexPos, m_forcestore))
    return false;
 
  if ((hasInnerHitStatus == 0) or m_forcestore)
    return true;

  bool failflag = false;
  const auto& trackHypotheses = getTrackHypotheses(v0Hypothesis);
  const int pdg_trackPlus   = trackPlus->getTrackFitResultWithClosestMass(
                                trackHypotheses.first)->getParticleType().getPDGCode();
  const int pdg_trackMinus  = trackMinus->getTrackFitResultWithClosestMass(
                                trackHypotheses.second)->getParticleType().getPDGCode();
 
  RecoTrack* recoTrackPlus_forRefit = nullptr;
  RecoTrack* recoTrackMinus_forRefit = nullptr;
  RecoTrack* cache_recoTrackPlus  = nullptr;
  RecoTrack* cache_recoTrackMinus = nullptr;
 
  unsigned int count_removeInnerHits = 0;
  while (hasInnerHitStatus != 0) {
    ++count_removeInnerHits;

    if (hasInnerHitStatus & 0x1) {
      // create a copy of the original RecoTrack w/o track fit
      recoTrackPlus_forRefit  = copyRecoTrack(recoTrackPlus);
      if (recoTrackPlus_forRefit == nullptr)
        return false;
      if (!cache_recoTrackPlus) {
        if (not removeInnerHits(recoTrackPlus, recoTrackPlus_forRefit, pdg_trackPlus, vertexPos)) {
          failflag = true;
          break;
        }
      } else {
        if (not removeInnerHits(cache_recoTrackPlus, recoTrackPlus_forRefit, pdg_trackPlus, vertexPos)) {
          failflag = true;
          break;
        }
      }
      cache_recoTrackPlus = recoTrackPlus_forRefit;
    } else if (recoTrackPlus_forRefit == nullptr) {
      // create a copy of the original RecoTrack w/ track fit (only once)
      recoTrackPlus_forRefit  = copyRecoTrackAndFit(recoTrackPlus, pdg_trackPlus);
      if (recoTrackPlus_forRefit == nullptr)
        return false;
    }

    if (hasInnerHitStatus & 0x2) {
      // create a copy of the original RecoTrack w/o track fit
      recoTrackMinus_forRefit = copyRecoTrack(recoTrackMinus);
      if (recoTrackMinus_forRefit == nullptr)
        return false;
      if (!cache_recoTrackMinus) {
        if (not removeInnerHits(recoTrackMinus, recoTrackMinus_forRefit, pdg_trackMinus, vertexPos)) {
          failflag = true;
          break;
        }
      } else {
        if (not removeInnerHits(cache_recoTrackMinus, recoTrackMinus_forRefit, pdg_trackMinus, vertexPos)) {
          failflag = true;
          break;
        }
      }
      cache_recoTrackMinus = recoTrackMinus_forRefit;
    } else if (recoTrackMinus_forRefit == nullptr) {
      // create a copy of the original RecoTrack w/ track fit (only once)
      recoTrackMinus_forRefit = copyRecoTrackAndFit(recoTrackMinus, pdg_trackMinus);
      if (recoTrackMinus_forRefit == nullptr)
        return false;
    }

    hasInnerHitStatus = 0;
    if (not vertexFitWithRecoTracks(trackPlus, trackMinus, recoTrackPlus_forRefit, recoTrackMinus_forRefit,
                                    v0Hypothesis, hasInnerHitStatus, vertexPos, false)) {
      B2DEBUG(22, "Vertex refit failed, or rejected by invariant mass cut.");
      failflag = true;
      break;
    } else if (hasInnerHitStatus == 0)
      isHitRemoved = true;
    if (count_removeInnerHits >= 5) {
      B2WARNING("Inner hits remained after " << count_removeInnerHits << " times of removing inner hits!");
      failflag = true;
      break;
    }
  }

  if (failflag) {
    bool forcestore = true;
    if (not vertexFitWithRecoTracks(trackPlus, trackMinus, recoTrackPlus, recoTrackMinus, v0Hypothesis,
                                    hasInnerHitStatus, vertexPos, forcestore)) {
      B2DEBUG(22, "Original vertex fit fails. Possibly rejected by invariant mass cut.");
      return false;
    } else
      isForceStored = true;
  }
 
  return true;
}

bool V0Fitter::vertexFitWithRecoTracks(const Track* trackPlus, const Track* trackMinus,
                                       RecoTrack* recoTrackPlus, RecoTrack* recoTrackMinus,
                                       const Const::ParticleType& v0Hypothesis,
                                       unsigned int& hasInnerHitStatus, ROOT::Math::XYZVector& vertexPos,
                                       const bool forceStore)
{
  const auto& trackHypotheses = getTrackHypotheses(v0Hypothesis);
 
  // make a clone, not use the reference so that the genfit::Track and its TrackReps will not be altered.
  genfit::Track gfTrackPlus = RecoTrackGenfitAccess::getGenfitTrack(*recoTrackPlus);
  const int pdgTrackPlus = trackPlus->getTrackFitResultWithClosestMass(trackHypotheses.first)->getParticleType().getPDGCode();
  const genfit::AbsTrackRep* plusRepresentation = recoTrackPlus->getTrackRepresentationForPDG(pdgTrackPlus);
  if ((plusRepresentation == nullptr) or (not recoTrackPlus->wasFitSuccessful(plusRepresentation))) {
    B2ERROR("Track hypothesis with closest mass not available. Should never happen, but I can continue safely anyway.");
    return false;
  }
 
  // make a clone, not use the reference so that the genfit::Track and its TrackReps will not be altered.
  genfit::Track gfTrackMinus = RecoTrackGenfitAccess::getGenfitTrack(*recoTrackMinus);
  const int pdgTrackMinus = trackMinus->getTrackFitResultWithClosestMass(trackHypotheses.second)->getParticleType().getPDGCode();
  const genfit::AbsTrackRep* minusRepresentation = recoTrackMinus->getTrackRepresentationForPDG(pdgTrackMinus);
  if ((minusRepresentation == nullptr) or (not recoTrackMinus->wasFitSuccessful(minusRepresentation))) {
    B2ERROR("Track hypothesis with closest mass not available. Should never happen, but I can continue safely anyway.");
    return false;
  }

  const std::vector<genfit::AbsTrackRep*>& repsPlus = gfTrackPlus.getTrackReps();
  const std::vector<genfit::AbsTrackRep*>& repsMinus = gfTrackMinus.getTrackReps();
  if (repsPlus.size() == repsMinus.size()) {
    for (unsigned int id = 0; id < repsPlus.size(); id++) {
      if (abs(repsPlus[id]->getPDG()) == pdgTrackPlus)
        gfTrackPlus.setCardinalRep(id);
      if (abs(repsMinus[id]->getPDG()) == pdgTrackMinus)
        gfTrackMinus.setCardinalRep(id);
    }
  }
  else {
    for (unsigned int id = 0; id < repsPlus.size(); id++) {
      if (abs(repsPlus[id]->getPDG()) == pdgTrackPlus)
        gfTrackPlus.setCardinalRep(id);
    }
    for (unsigned int id = 0; id < repsMinus.size(); id++) {
      if (abs(repsMinus[id]->getPDG()) == pdgTrackMinus)
        gfTrackMinus.setCardinalRep(id);
    }
  }

  genfit::MeasuredStateOnPlane stPlus  = recoTrackPlus->getMeasuredStateOnPlaneFromFirstHit(plusRepresentation);
  genfit::MeasuredStateOnPlane stMinus = recoTrackMinus->getMeasuredStateOnPlaneFromFirstHit(minusRepresentation);
 
  genfit::GFRaveVertex vert;
  if (not fitGFRaveVertex(gfTrackPlus, gfTrackMinus, vert)) {
    return false;
  }
 
  const ROOT::Math::XYZVector& posVert = ROOT::Math::XYZVector(vert.getPos());
  vertexPos = posVert;

  if (posVert.Rho() < m_beamPipeRadius) {
    return false;
  } else {
    if (vert.getChi2() > m_vertexChi2CutOutside) {
      B2DEBUG(22, "Vertex outside beam pipe, chi^2 too large.");
      return false;
    }
  }
 
  B2DEBUG(22, "Vertex accepted.");
 
  if (not extrapolateToVertex(stPlus, stMinus, posVert, hasInnerHitStatus)) {
    return false;
  }

  if (forceStore || hasInnerHitStatus == 0) {
    const genfit::GFRaveTrackParameters* tr0 = vert.getParameters(0);
    const genfit::GFRaveTrackParameters* tr1 = vert.getParameters(1);
    ROOT::Math::PxPyPzMVector lv0(tr0->getMom().Px(), tr0->getMom().Py(), tr0->getMom().Pz(), trackHypotheses.first.getMass());
    ROOT::Math::PxPyPzMVector lv1(tr1->getMom().Px(), tr1->getMom().Py(), tr1->getMom().Pz(), trackHypotheses.second.getMass());
    double v0InvMass = (lv0 + lv1).M();
    if (v0Hypothesis == Const::Kshort) {
      if (v0InvMass < std::get<0>(m_invMassRangeKshort) || v0InvMass > std::get<1>(m_invMassRangeKshort)) {
        B2DEBUG(22, "Kshort vertex rejected, invariant mass out of range.");
        return false;
      }
    } else if (v0Hypothesis == Const::Lambda || v0Hypothesis == Const::antiLambda) {
      if (v0InvMass < std::get<0>(m_invMassRangeLambda) || v0InvMass > std::get<1>(m_invMassRangeLambda)) {
        B2DEBUG(22, "Lambda vertex rejected, invariant mass out of range.");
        return false;
      }
    } else if (v0Hypothesis == Const::photon) {
      if (v0InvMass < std::get<0>(m_invMassRangePhoton) || v0InvMass > std::get<1>(m_invMassRangePhoton)) {
        B2DEBUG(22, "Photon vertex rejected, invariant mass out of range.");
        return false;
      }
    }
 

    const double Bz = BFieldManager::getFieldInTesla({0, 0, 0}).Z();
 
    int sharedInnermostCluster = checkSharedInnermostCluster(recoTrackPlus, recoTrackMinus);
 
    TrackFitResult* tfrPlusVtx =  buildTrackFitResult(gfTrackPlus, recoTrackPlus, stPlus,  Bz, trackHypotheses.first,
                                                      sharedInnermostCluster);
    TrackFitResult* tfrMinusVtx = buildTrackFitResult(gfTrackMinus, recoTrackMinus, stMinus, Bz, trackHypotheses.second,
                                                      sharedInnermostCluster);
 
    B2DEBUG(20, "Creating new V0.");
    const auto* v0 = m_v0s.appendNew(std::make_pair(trackPlus, tfrPlusVtx),
                                     std::make_pair(trackMinus, tfrMinusVtx),
                                     posVert.X(), posVert.Y(), posVert.Z());
 
    if (m_validation) {
      B2DEBUG(24, "Create StoreArray and Output for validation.");
      const auto* validationV0 = m_validationV0s.appendNew(
                                   std::make_pair(trackPlus, tfrPlusVtx),
                                   std::make_pair(trackMinus, tfrMinusVtx),
                                   ROOT::Math::XYZVector(vert.getPos()),
                                   vert.getCov(),
                                   (lv0 + lv1).P(),
                                   (lv0 + lv1).M(),
                                   vert.getChi2()
                                 );
      v0->addRelationTo(validationV0);
    }
  }
  return true;
}
 
RecoTrack* V0Fitter::copyRecoTrack(const RecoTrack* origRecoTrack)
{
  RecoTrack* newRecoTrack = origRecoTrack->copyToStoreArray(m_copiedRecoTracks);
  newRecoTrack->addHitsFromRecoTrack(origRecoTrack);
  newRecoTrack->addRelationTo(origRecoTrack);
  return newRecoTrack;
}
 
RecoTrack* V0Fitter::copyRecoTrackAndFit(const RecoTrack* origRecoTrack, const int trackPDG)
{
  Const::ChargedStable particleUsedForFitting(std::abs(trackPDG));
  const genfit::AbsTrackRep* origTrackRep = origRecoTrack->getTrackRepresentationForPDG(std::abs(
                                              trackPDG));
 
  RecoTrack* newRecoTrack = copyRecoTrack(origRecoTrack);

  TrackFitter fitter;
  if (not fitter.fit(*newRecoTrack, particleUsedForFitting)) {
    // This is not expected, but happens sometimes.
    B2DEBUG(20, "track fit failed for copied RecoTrack.");
    if (not origRecoTrack->wasFitSuccessful(origTrackRep))
      B2DEBUG(20, "\t original track fit was also failed.");
    return nullptr;
  }
 
  return newRecoTrack;
}
 
bool V0Fitter::removeInnerHits(RecoTrack* prevRecoTrack, RecoTrack* recoTrack,
                               const int trackPDG, const ROOT::Math::XYZVector& vertexPosition)
{
  if (!prevRecoTrack || !recoTrack) {
    B2ERROR("Input recotrack is nullptr!");
    return false;
  }
  Const::ChargedStable particleUsedForFitting(std::abs(trackPDG));
  const genfit::AbsTrackRep* prevTrackRep = prevRecoTrack->getTrackRepresentationForPDG(std::abs(
                                              trackPDG));

  const std::vector<RecoHitInformation*>& recoHitInformations = recoTrack->getRecoHitInformations(true);
  const std::vector<RecoHitInformation*>& prevRecoHitInformations = prevRecoTrack->getRecoHitInformations(
        true);
  unsigned int nRemoveHits = 0;
  if (recoHitInformations.size() != prevRecoHitInformations.size()) {
    B2WARNING("Copied RecoTrack has different number of hits from its original RecoTrack!");
    return false;
  }
  for (nRemoveHits = 0; nRemoveHits < recoHitInformations.size(); ++nRemoveHits) {
    if (!prevRecoHitInformations[nRemoveHits]->useInFit()) {
      recoHitInformations[nRemoveHits]->setUseInFit(false);
      continue;
    }
    try {
      genfit::MeasuredStateOnPlane stPrevRecoHit = prevRecoTrack->getMeasuredStateOnPlaneFromRecoHit(
                                                     prevRecoHitInformations[nRemoveHits]);
      double extralength = stPrevRecoHit.extrapolateToPoint(XYZToTVector(vertexPosition));
      if (extralength > 0) {
        recoHitInformations[nRemoveHits]->setUseInFit(false);
      } else
        break;
    } catch (NoTrackFitResult()) {
      B2WARNING("Exception: no FitterInfo assigned for TrackPoint created from this RecoHit.");
      recoHitInformations[nRemoveHits]->setUseInFit(false);
      continue;
    } catch (...) {
      B2DEBUG(22, "Could not extrapolate track to vertex when removing inner hits, aborting.");
      return false;
    }
  }
 
  if (nRemoveHits == 0) {
    // This is not expected, but can happen if the track fit is different between copied and original RecoTrack.
    B2DEBUG(20, "No hits removed in removeInnerHits, aborted. Switching to use the original RecoTrack.");
    return false;
  }
 
  if (recoHitInformations.size() <= nRemoveHits) {
    B2DEBUG(20, "Removed all the RecoHits in the RecoTrack, aborted. Switching to use the original RecoTrack.");
    return false;
  }

  if (recoHitInformations[nRemoveHits - 1]->getTrackingDetector() == RecoHitInformation::RecoHitDetector::c_SVD) {
    if (recoHitInformations[nRemoveHits]->getTrackingDetector() == RecoHitInformation::RecoHitDetector::c_SVD) {
      const SVDCluster* lastRemovedSVDHit = recoHitInformations[nRemoveHits - 1]->getRelatedTo<SVDCluster>();
      const SVDCluster* nextSVDHit        = recoHitInformations[nRemoveHits]  ->getRelatedTo<SVDCluster>();
      if (!lastRemovedSVDHit || !nextSVDHit) B2ERROR("Last/Next SVD hit is null!");
      else {
        if (lastRemovedSVDHit->getSensorID() == nextSVDHit->getSensorID() &&
            lastRemovedSVDHit->isUCluster() && !(nextSVDHit->isUCluster())) {
          recoHitInformations[nRemoveHits]->setUseInFit(false);
          ++nRemoveHits;
        }
      }
    }
  }

  if (!m_useOnlyOneSVDHitPair &&
      recoHitInformations[nRemoveHits - 1]->getTrackingDetector() == RecoHitInformation::RecoHitDetector::c_SVD &&
      recoHitInformations.size() > nRemoveHits + 2) {
    if (recoHitInformations[nRemoveHits  ]->getTrackingDetector() == RecoHitInformation::RecoHitDetector::c_SVD && 
        recoHitInformations[nRemoveHits + 1]->getTrackingDetector() == RecoHitInformation::RecoHitDetector::c_SVD && 
        recoHitInformations[nRemoveHits + 2]->getTrackingDetector() != RecoHitInformation::RecoHitDetector::c_SVD) { 
      recoHitInformations[nRemoveHits  ]->setUseInFit(false);
      recoHitInformations[nRemoveHits + 1]->setUseInFit(false);
      nRemoveHits += 2;
    }
  }
 
  B2DEBUG(22, nRemoveHits << " inner hits removed.");

  TrackFitter fitter;
  if (not fitter.fit(*recoTrack, particleUsedForFitting)) {
    B2DEBUG(20, "track fit failed after removing inner hits.");
    if (not prevRecoTrack->wasFitSuccessful(prevTrackRep))
      B2DEBUG(20, "\t previous track fit was also failed.");
    return false;
  }
 
  return true;
}
 
int V0Fitter::checkSharedInnermostCluster(const RecoTrack* recoTrackPlus, const RecoTrack* recoTrackMinus)
{
  int flag = 0; // -1 for exception
 
  // get the innermost hit for plus/minus-daughter
  const std::vector<RecoHitInformation*>& recoHitInformationsPlus  = recoTrackPlus->getRecoHitInformations(
        true); // true for sorted info.
  const std::vector<RecoHitInformation*>& recoHitInformationsMinus = recoTrackMinus->getRecoHitInformations(
        true);// true for sorted info.
  unsigned int iInnermostHitPlus, iInnermostHitMinus;
  for (iInnermostHitPlus = 0 ; iInnermostHitPlus < recoHitInformationsPlus.size() ; ++iInnermostHitPlus)
    if (recoHitInformationsPlus[iInnermostHitPlus]->useInFit()) break;
  for (iInnermostHitMinus = 0 ; iInnermostHitMinus < recoHitInformationsMinus.size() ; ++iInnermostHitMinus)
    if (recoHitInformationsMinus[iInnermostHitMinus]->useInFit()) break;
  if (iInnermostHitPlus == recoHitInformationsPlus.size() || iInnermostHitMinus == recoHitInformationsMinus.size()) {
    B2WARNING("checkSharedInnermostCluster function called for recoTrack including no hit used for fit! This should not happen!");
    return -1;
  }
  const auto& recoHitInfoPlus  = recoHitInformationsPlus[iInnermostHitPlus];
  const auto& recoHitInfoMinus = recoHitInformationsMinus[iInnermostHitMinus];
 
  if (recoHitInfoPlus->getTrackingDetector() == recoHitInfoMinus->getTrackingDetector()) {
    if (recoHitInfoPlus->getTrackingDetector() == RecoHitInformation::c_PXD) {
      const PXDCluster* clusterPlus = recoHitInfoPlus->getRelatedTo<PXDCluster>();
      const PXDCluster* clusterMinus = recoHitInfoMinus->getRelatedTo<PXDCluster>();
      if (clusterPlus == clusterMinus) { // if they share a same PXDCluster, set the flag
        flag = 3; // PXD cluster is a 2D-hit
      }
    } else if (recoHitInfoPlus->getTrackingDetector() == RecoHitInformation::c_SVD) {
      const SVDCluster* clusterPlus = recoHitInfoPlus->getRelatedTo<SVDCluster>();
      const SVDCluster* clusterMinus = recoHitInfoMinus->getRelatedTo<SVDCluster>();
      if (clusterPlus->isUCluster() && clusterMinus->isUCluster()) {
        if (recoHitInformationsPlus.size() > iInnermostHitPlus + 1
            && recoHitInformationsMinus.size() > iInnermostHitMinus + 1) { // not to access an array out of boundary
          const auto& recoHitInfoNextPlus  = recoHitInformationsPlus[iInnermostHitPlus + 1];
          const auto& recoHitInfoNextMinus = recoHitInformationsMinus[iInnermostHitMinus + 1];
          // sanity check to access next hits
          if (recoHitInfoNextPlus->useInFit() && recoHitInfoNextMinus->useInFit()  // this should be satisfied by default
              && recoHitInfoNextPlus->getTrackingDetector()  == RecoHitInformation::c_SVD
              && recoHitInfoNextMinus->getTrackingDetector() == RecoHitInformation::c_SVD) {
            const SVDCluster* clusterNextPlus = recoHitInfoNextPlus->getRelatedTo<SVDCluster>();
            const SVDCluster* clusterNextMinus = recoHitInfoNextMinus->getRelatedTo<SVDCluster>();
            if (!(clusterNextPlus->isUCluster()) && !(clusterNextMinus->isUCluster())
                && clusterPlus->getSensorID() == clusterNextPlus->getSensorID()
                && clusterMinus->getSensorID() == clusterNextMinus->getSensorID()) {
              if (clusterPlus == clusterMinus)
                flag += 2; // SVD U-cluster is shared
              if (clusterNextPlus == clusterNextMinus)
                flag += 1; // SVD V-cluster is shared
            } else {
              B2WARNING("SVD cluster to be paired is not on V-side, or not on the same sensor.");
              return -1;
            }
          } else {
            B2WARNING("No SVD cluster to be paired.");
            return -1;
          }
        } else {
          B2WARNING("Innermost SVD U-cluster is the only hit in a daughter track. This should not happen.");
          return -1;
        }
      } else {
        B2WARNING("No SVD U-cluster in the innermost cluster.");
        return -1;
      }
    }
  }
  return flag;
 
}