SiTrajectory_xk.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
#include "SiSPSeededTrackFinderData/SiTrajectory_xk.h"
 
#include <format>
#include <iostream>

// Set work information to trajectory
 
void InDet::SiTrajectory_xk::setTools(const InDet::SiTools_xk* t)
{
  m_tools = t;
  for (int i=0; i!=300; ++i) m_elements[i].setTools(t);
} 
 
void InDet::SiTrajectory_xk::setParameters()
{
  for (int i=0; i!=300; ++i) m_elements[i].setParameters();
} 

// Erase trajector element
 
void InDet::SiTrajectory_xk::erase(int n) 
{
  if (n>=0 && n<m_nElements) {
    for (int i=n; i!=m_nElements-1; ++i) m_elementsMap[i] = m_elementsMap[i+1];
    --m_nElements;
  }
}

// Trajectory conversion to TrackStateOnSurface  
 
Trk::TrackStates
InDet::SiTrajectory_xk::convertToTrackStateOnSurface(int cosmic, const EventContext& ctx)
{
  if (!cosmic ||  m_elements[m_elementsMap[m_firstElement]].parametersUB().parameters()[2] < 0.) {
    return convertToTrackStateOnSurface(ctx);
  }
  return convertToTrackStateOnSurfaceWithNewDirection(ctx);
}

// Trajectory conversion to TrackStateOnSurface  with old direction
 
Trk::TrackStates
InDet::SiTrajectory_xk::convertToTrackStateOnSurface(const EventContext& ctx)
{
 
  auto dtsos = Trk::TrackStates();
 
  bool multi = m_tools->multiTrack();
  if (m_nclusters <= m_tools->clustersmin() ||
      pTfirst() < m_tools->pTmin()) multi = false;
 
  int i = m_firstElement;
  
  const Trk::TrackStateOnSurface* 
    tsos = m_elements[m_elementsMap[i]].trackStateOnSurface(false,true,multi,1,ctx);
 
  if (tsos) dtsos.push_back(tsos);
 
  for (++i; i!=m_lastElement; ++i) {
 
    int m = m_elementsMap[i];
    if (m_elements[m].cluster() || m_elements[m].clusterNoAdd() ) {
      tsos = m_elements[m].trackStateOnSurface(false,false,multi,0,ctx);
      if (tsos) dtsos.push_back(tsos);
    }
  }
 
  i = m_lastElement;
  tsos = m_elements[m_elementsMap[i]].trackStateOnSurface(false,false,multi,2,ctx);
  if (tsos) dtsos.push_back(tsos);
 
  if (multi) {
    m_ntos = 0;
    for (int i=m_firstElement; i<=m_lastElement; ++i) {
      
      int m = m_elementsMap[i];
      if (!m_elements[m].ntsos()) continue;
      m_atos[m_ntos  ] = m;
      m_itos[m_ntos++] = 0;
    }
  }
  return dtsos;
}

// Trajectory conversion to TrackStateOnSurface  with new direction
 
Trk::TrackStates
InDet::SiTrajectory_xk::convertToTrackStateOnSurfaceWithNewDirection(const EventContext& ctx)
{
 
  auto dtsos = Trk::TrackStates();
 
  bool multi = m_tools->multiTrack();
  if (pTfirst() < m_tools->pTmin()) multi = false;
 
  int i = m_lastElement;
 
  const Trk::TrackStateOnSurface* 
    tsos = m_elements[m_elementsMap[i]].trackStateOnSurface(true,true,multi,2,ctx);
 
  if (tsos) dtsos.push_back(tsos);
 
  for (--i; i!=m_firstElement; --i) {
 
    int m = m_elementsMap[i];
    if (m_elements[m].cluster() || m_elements[m].clusterNoAdd() ) {
      tsos = m_elements[m].trackStateOnSurface(true,false,multi,0,ctx);
      if (tsos) dtsos.push_back(tsos);
    }
  }
 
  i = m_firstElement;
  tsos = m_elements[m_elementsMap[i]].trackStateOnSurface(true,false,multi,1,ctx);
  if (tsos) dtsos.push_back(tsos);
 
  return dtsos;
}

// Trajectory conversion to simple TrackStateOnSurface   
 
Trk::TrackStates
InDet::SiTrajectory_xk::convertToSimpleTrackStateOnSurface(int cosmic, const EventContext& ctx)
{
  if (!cosmic ||  m_elements[m_elementsMap[m_firstElement]].parametersUB().parameters()[2] < 0.) {
    return convertToSimpleTrackStateOnSurface(ctx);
  }
  return convertToSimpleTrackStateOnSurfaceWithNewDirection();
}

// Trajectory conversion to simple TrackStateOnSurface  with old direction
 
Trk::TrackStates 
InDet::SiTrajectory_xk::convertToSimpleTrackStateOnSurface(const EventContext& ctx)
{
  auto dtsos = Trk::TrackStates();
 
  int i = m_firstElement;
  
  const Trk::TrackStateOnSurface* 
    tsos = m_elements[m_elementsMap[i]].trackPerigeeStateOnSurface(ctx);
 
  if (tsos) dtsos.push_back(tsos);
  
  tsos = m_elements[m_elementsMap[i]].trackSimpleStateOnSurface(false,m_tools->useFastTracking(),m_tools->useFastTracking());
 
  if (tsos) dtsos.push_back(tsos);
  
  int lastClusterElement = 0;
  for (int j=m_lastElement; j>=i; j--) {
     int m = m_elementsMap[j];
     if (m_elements[m].cluster()) { 
    lastClusterElement = j;
    break;
     }
  }
  if( lastClusterElement==0 || lastClusterElement==i ) return dtsos;
 
  for (++i; i<std::min(lastClusterElement,m_lastElement); ++i) {
    
    int m = m_elementsMap[i];
    if (m_elements[m].cluster()) {
      tsos = m_elements[m].trackSimpleStateOnSurface(false,false,0);
      if (tsos) dtsos.push_back(tsos);
    }
  }
 
  i = std::min(lastClusterElement,m_lastElement);
  tsos = m_elements[m_elementsMap[i]].trackSimpleStateOnSurface(false,true,2);
  if (tsos) dtsos.push_back(tsos);
 
  return dtsos;
}

// Trajectory conversion to simple TrackStateOnSurface with new direction
 
Trk::TrackStates
InDet::SiTrajectory_xk::convertToSimpleTrackStateOnSurfaceWithNewDirection()
{
  auto dtsos = Trk::TrackStates();
 
  int i = m_lastElement;
 
  const Trk::TrackStateOnSurface* 
    tsos = m_elements[m_elementsMap[i]].trackSimpleStateOnSurface(true,true,2);
 
  if (tsos) dtsos.push_back(tsos);
 
  for (--i; i!=m_firstElement; --i) {
 
    int m = m_elementsMap[i];
    if (m_elements[m].cluster() || m_elements[m].clusterNoAdd() ) {
      tsos = m_elements[m].trackSimpleStateOnSurface(true,false,0);
      if (tsos) dtsos.push_back(tsos);
    }
  }
 
  i = m_firstElement;
  tsos = m_elements[m_elementsMap[i]].trackSimpleStateOnSurface(true,false,1);
  if (tsos) dtsos.push_back(tsos);
 
  return dtsos;
}

// Trajectory conversion to simple TrackStateOnSurface   
// Only for Disappearing Track Trigger that uses also failed tracks
 
Trk::TrackStates
InDet::SiTrajectory_xk::convertToSimpleTrackStateOnSurfaceForDisTrackTrigger(int cosmic, const EventContext& ctx)
{
  if (!cosmic ||  m_elements[m_elementsMap[m_firstElement]].parametersUB().parameters()[2] < 0.) {
    return convertToSimpleTrackStateOnSurfaceForDisTrackTrigger(ctx);
  }
  return convertToSimpleTrackStateOnSurfaceWithNewDirection();
}

// Trajectory conversion to simple TrackStateOnSurface  with old direction
// Only for Disappearing Track Trigger that uses also failed tracks
 
Trk::TrackStates 
InDet::SiTrajectory_xk::convertToSimpleTrackStateOnSurfaceForDisTrackTrigger(const EventContext& ctx)
{
  auto dtsos = Trk::TrackStates();
 
  int i = m_firstElement;
  
  const Trk::TrackStateOnSurface* 
    tsos = m_elements[m_elementsMap[i]].trackPerigeeStateOnSurface(ctx);
 
  if (tsos) dtsos.push_back(tsos);
  
  tsos = m_elements[m_elementsMap[i]].trackSimpleStateOnSurface(false,false,0);
 
  if (tsos) dtsos.push_back(tsos);
  
  int lastClusterElement = 0;
  for (int j=m_lastElement; j>=i; j--) {
     int m = m_elementsMap[j];
     if (m_elements[m].cluster()) { 
    lastClusterElement = j;
    break;
     }
  }
  if( lastClusterElement==0 || lastClusterElement==i ) return dtsos;
 
  for (++i; i<std::min(lastClusterElement,m_lastElement); ++i) {
    
    int m = m_elementsMap[i];
    if (m_elements[m].cluster()) {
      tsos = m_elements[m].trackSimpleStateOnSurface(false,false,0);
      if (tsos) dtsos.push_back(tsos);
    }
  }
 
  i = std::min(lastClusterElement,m_lastElement);
  tsos = m_elements[m_elementsMap[i]].trackSimpleStateOnSurface(false,true,2);
  if (tsos) dtsos.push_back(tsos);
 
  return dtsos;
}

// FitQuality production
 
std::unique_ptr<Trk::FitQuality> InDet::SiTrajectory_xk::convertToFitQuality() const{
  double xi2 = m_elements[m_elementsMap[m_firstElement]].xi2totalB();
  return std::make_unique<Trk::FitQuality>(xi2, (m_ndf - 5));
}

// Test is it new track
 
bool InDet::SiTrajectory_xk::isNewTrack
(std::multimap<const Trk::PrepRawData*,const Trk::Track*>& map) const
{
  if (m_firstElement==-100) return false;//i.e. the int array never had elements inserted
  const Trk::PrepRawData* prd   [100];
  std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator 
    ti,t[100],te = map.end();
  int n = 0 ;
  for (int i=m_firstElement; i<=m_lastElement; ++i) {
    if (n >= 100) break;
    int m = m_elementsMap[i];
    
    if (m_elements[m].cluster()) {
      prd[n] = m_elements[m].cluster();
      t[n] = map.find(prd[n]);
      if (t[n]==te) return true;
      ++n;
    } else if (m_elements[m].clusterNoAdd()) {
      prd[n] = m_elements[m].clusterNoAdd();
      t[n] = map.find(prd[n]);
      if (t[n]==te) return true;
      ++n;
    }
  }
 
  int nclt = m_nclusters + m_nclustersNoAdd;
  
  for (int i=0; i!=n; ++i) {
    int nclmax = 0;
    for (ti=t[i]; ti!=te; ++ti) {
      if ( (*ti).first != prd[i] ) break;
      int ncl = (*ti).second->measurementsOnTrack()->size();
      if (ncl > nclmax) nclmax = ncl;
    }   
    if (nclt > nclmax) return true;
  }
  return false;
}

// Overload of << operator std::ostream
 
std::ostream& InDet::operator << 
(std::ostream& sl,const InDet::SiTrajectory_xk& se)
{ 
  return se.dump(sl);
}   

// Dumps relevant information into the ostream
 
std::ostream& InDet::SiTrajectory_xk::dump( std::ostream& out ) const
{
  if (m_nElements <=0 ) {
    out<<"Trajectory does not exist"<<std::endl;
    return out;
  }
  if (m_firstElement >= m_lastElement ) {
    out<<"Trajectory is wrong"<<std::endl;
    return out;
  }
 
  out<<"|--------------------------------------------------------------------------------------------------------|"
     <<std::endl;
  out<<"|                       TRAJECTORY "
     <<"                                                                      |"
     <<std::endl;
 
    out << std::format("| Has{:>3} ({:>3}) elements and {:>2} ({:>2}) clusters and {:>2} weighted clusters and quality = {:>12.5g}         |\n",
                m_nElements,
                m_nActiveElements,
                m_nclusters + m_nclustersNoAdd,
                m_nclustersNoAdd,
                m_ndf,
                quality());
    out << std::format("| Has number of holes before, inside, after and gap= {:>2}{:>2}{:>2}{:>3}                                           |\n",
                m_nHolesBefore,
                m_nholes,
                m_nHolesAfter,
                m_dholes);
  out<<"|                                                                                        F   B           |"
     <<std::endl;
 
  out<<"|---|--|---|-----|-|-|---------|---------|---------|---------|----------|---------|-|--|--|--|-|-|-|-|-|-|---------|"
     <<std::endl;
  out<<"| # |SS|  D| Ncl |C|O|   Xi2F  |   Xi2B  | Az.Angle| Radius  |  pT(GeV) |  dZ/dR  |N|In|Lf|Lb|S|D|H|G|H|G|   Step  |"
     <<std::endl;
  out<<"|---|--|---|-----|-|-|---------|---------|---------|---------|----------|---------|-|--|--|--|-|-|-|-|-|-|---------|"
     <<std::endl;
 
  for (int i=0; i!=m_nElements; ++i) {
    
    int m = m_elementsMap[i];
 
    std::string DET = "D ";
    const InDetDD::SiDetectorElement* D = m_elements[m].detElement();
    std::string DE = " ";
    if (m_elements[m].detstatus() < 0) DE = "-";
    if (D) {
      if (D->isPixel()) {
        if (D->isBarrel()) DET = "Pb"; else DET = "Pe";
      }
      else if (D->isSCT()) {
        if (D->isBarrel()) DET = "Sb"; else DET = "Se";
      }
    }
    int c = 0;
    if (m_elements[m].detstatus() > 0) c = m_elements[m].numberClusters();
 
    out << std::format("|{:>3}", static_cast<unsigned>(i));
 
    std::string S0="  ";
    if (m_firstElement == i) S0="=>";
    if (m_lastElement  == i) S0="=>";
 
    std::string S1=" ";
    std::string S2=" ";
    if (m_elements[m].cluster     ()) S1="+";
    if (m_elements[m].clusterNoAdd()) S2="+";
 
    out << std::format("|{}|{:>1}{:>2}|{:>5}|{}|{}|", S0, DE, DET, c, S1, S2);
 
    if (m_elements[m].status()) {
 
      out << std::format("{:>9.3g}|{:>9.3g}|", m_elements[m].xi2F(), m_elements[m].xi2B());
 
      double ra = 0.;
      double pt = 0.;
      double tz = 0.;
      double fa = 0.;
 
      if (m_elements[m].status()==1) {
 
        if (m_elements[m].cluster()) {
 
          Amg::Vector3D gp = m_elements[m].parametersUF().position();
          ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
          fa = atan2(gp.y(),gp.x());
          pt = m_elements[m].parametersUF().momentum().perp();
          tz = m_elements[m].parametersUF().cotTheta();
        }
        else {
 
          Amg::Vector3D gp = m_elements[m].parametersPF().position();
          ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
          fa = atan2(gp.y(),gp.x());
          pt = m_elements[m].parametersPF().momentum().perp();
          tz = m_elements[m].parametersPF().cotTheta();
        }
      }
      else if ((m_tools->useFastTracking() and m_elements[m].status()>2) or m_elements[m].status()==2) {
 
        if (m_elements[m].cluster()) {
 
          Amg::Vector3D gp = m_elements[m].parametersUB().position();
          ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
          fa = atan2(gp.y(),gp.x());
          pt = m_elements[m].parametersUB().momentum().perp();
          tz = m_elements[m].parametersUB().cotTheta();
        }
        else {
 
          Amg::Vector3D gp = m_elements[m].parametersPB().position();
          ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
          fa = atan2(gp.y(),gp.x());
          pt = m_elements[m].parametersPB().momentum().perp();
          tz = m_elements[m].parametersPB().cotTheta();
        }
      }
      else {
 
        Trk::PatternTrackParameters S1,SM,S2(m_elements[m].parametersPF());
 
        if (m_elements[m].cluster()) S1 = m_elements[m].parametersUB();
        else                        S1 = m_elements[m].parametersPB();
 
        bool QA = m_tools->updatorTool()->combineStates(S1,S2,SM);
 
        if (QA) {
 
          Amg::Vector3D gp = SM.position();
          ra = sqrt(gp.x()*gp.x()+gp.y()*gp.y());
          fa = atan2(gp.y(),gp.x());
          pt = SM.momentum().perp();
          tz = SM.cotTheta();
        }
      }
      out << std::format("{:>9.4g}|{:>9.4g}|{:>10.4g}|{:>9.4g}|{:>1}|{:>2}|{:>2}|{:>2}|{:>1}|{:>1}|{:>1}|{:>1}|{:>1}|{:>1}|{:>9.4g}|",
             fa,
             ra,
             pt*0.001,
             tz,
             static_cast<unsigned>(m_elements[m].noiseModel()),
             m_elements[m].inside(),
             static_cast<unsigned>(m_elements[m].nlinksF()),
             static_cast<unsigned>(m_elements[m].nlinksB()),
             static_cast<unsigned>(m_elements[m].status()),
             static_cast<unsigned>(m_elements[m].difference()),
             static_cast<unsigned>(m_elements[m].nholesF()),
             static_cast<unsigned>(m_elements[m].dholesF()),
             static_cast<unsigned>(m_elements[m].nholesB()),
             static_cast<unsigned>(m_elements[m].dholesB()),
             m_elements[m].step());
    }
    else                       {
      out<<"         |";
      out<<"         |";
      out<<"         |";
      out<<"         |";
      out<<"          |";
      out<<"         |";
      out<<" |";
      out<<"  |";
      out<<"  |";
      out<<"  |";
      out<<" |";
      out<<" |";
      out<<" |";
      out<<" |";
      out<<" |";
      out<<" |";
      out << std::format("{:>9.4g}|", m_elements[m].step());
    }
    out<<std::endl;
  }
  out<<"|---|--|---|-----|-|-|---------|---------|---------|---------|----------|---------|-|--|--|--|-|-|-|-|-|-|---------|"
     <<std::endl;
  return out;
}   

// pT seed estimation
 
double InDet::SiTrajectory_xk::pTseed
 (const Trk::TrackParameters                            & Tp,
  std::vector<const InDet::SiCluster*>                    & Cl,
  std::vector<const InDet::SiDetElementBoundaryLink_xk*>& DE,
  const EventContext                                    & ctx)
{
  double Xi2cut     =  30.;
 
  InDet::SiClusterCollection::const_iterator  sib,sie;
  std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator r=DE.begin(),re=DE.end();
  std::vector<const InDet::SiCluster*>                    ::iterator s=Cl.begin();
 
  int n = 0;
  if(!m_elements[n].set(1,(*r),sib,sie,(*s),ctx) ) return 0.;
  if(!m_elements[n].firstTrajectorElement(Tp,ctx)) return 0.;
 
  for(++r; r!=re; ++r) {
    ++n; ++s;
    if(!m_elements[n].set(1,(*r),sib,sie,(*s),ctx)                    ) return 0.;
    if(!m_elements[n].ForwardPropagationWithoutSearch(m_elements[n-1], ctx)) return 0.;
    if( m_elements[n].xi2F()      >      Xi2cut                       ) return 0.;
  }
  return m_elements[n].parametersUF().momentum().perp();
}
 

// Initiate trajectory
 
bool InDet::SiTrajectory_xk::initialize
(bool PIX,
 bool SCT,
 const InDet::PixelClusterContainer*                  PIXc       ,
 const InDet::SCT_ClusterContainer*                   SCTc       ,
 const Trk::TrackParameters                          & Tp        ,
 std::vector<const InDet::SiCluster*>                  & lSiCluster, 
 std::vector<const InDet::SiDetElementBoundaryLink_xk*>& DE      ,
 bool                                                & rquality  ,
 const EventContext                                  & ctx       )
{
  m_nholes          =    0;
  m_nHolesBefore         =    0;
  m_nHolesAfter         =    0;
  m_dholes          =    0;
  m_nclusters       =    0;
  m_nclustersNoAdd  =    0;
  m_nElements       =    0;
  m_nActiveElements      =    0;
  m_firstElement    = -100;
  m_lastElement     =    0;
  m_ndfcut          =    0;
  rquality          = true;
  m_ntos            =    0;
  int    ndfwrong   =    0;
  double Xi2cut     = 2.*m_tools->xi2max();
 
  // radius of the dead cylinder
  double Rdead      = 142.5;
  // boolean to decide if initialisation is needed or not
  // initDeadMaterial is False (which means dead material needs be initialised)
  // for ITk fast tracking configuration
  bool initDeadMaterial = not(m_tools->isITkGeometry() and m_tools->useFastTracking());
 
  if(!initDeadMaterial and !m_surfacedead) m_surfacedead = std::make_unique<const Trk::CylinderSurface>(Rdead,5000.);
 
  std::vector<const InDet::SiCluster*>::iterator iter_cluster;
  if (lSiCluster.size() < 2) return false;
 
  std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator iter_boundaryLink,endBoundaryLinks=DE.end();
 
  int up    = 0;
  int last  = 0;

  for (iter_boundaryLink=DE.begin(); iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {

    const InDetDD::SiDetectorElement* detectorElement = (*iter_boundaryLink)->detElement();
    IdentifierHash           id = detectorElement->identifyHash();

    const InDet::SiCluster* theCluster = nullptr;

    if (detectorElement->isPixel()) {
      if (PIX) {  
 
    // Set dead material
    //
    // if already initialised, not doing it again
        if(not initDeadMaterial) {
          const Trk::PlaneSurface* pla = static_cast<const Trk::PlaneSurface*>(&detectorElement->surface());
          double R = pla->center().perp();
          if(R > Rdead) {
            initDeadMaterial = true;
            if(!m_elements[m_nElements].setDead(m_surfacedead.get())) return false;
            m_elementsMap[m_nElements] = m_nElements;
            if(m_nclusters && !lSiCluster.empty()) {
              if(!m_elements[m_nElements].ForwardPropagationWithoutSearch(m_elements[up], ctx)) return false;
              up = m_nElements;
            }
            if(++m_nElements==300) break;
          }
        }
 
        InDet::PixelClusterCollection::const_iterator iter_PixelClusterColl, iter_PixelClusterCollEnd;
        const InDet::PixelClusterCollection *clustersOnElement = (*PIXc).indexFindPtr(id);
        if (clustersOnElement!=nullptr && clustersOnElement->begin()!=clustersOnElement->end()) {

          iter_PixelClusterColl = clustersOnElement->begin();
          iter_PixelClusterCollEnd = clustersOnElement->end();
          
          for (iter_cluster=lSiCluster.begin(); iter_cluster!=lSiCluster.end(); ++iter_cluster) {
            if ((*iter_cluster)->detectorElement()==detectorElement) {
              if (m_nclusters==0){ 
                m_firstElement = m_nElements;
              }
              else{
                m_lastElement  = m_nElements;
              }
              ++m_nclusters;
              m_ndfcut+=2;
              theCluster=(*iter_cluster);
              iter_cluster=lSiCluster.erase(iter_cluster);
              break;
            }
          }
          bool valid_set = m_elements[m_nElements].set(1,(*iter_boundaryLink),iter_PixelClusterColl,iter_PixelClusterCollEnd,theCluster,ctx);
          if(m_tools->isITkGeometry() && !valid_set) return false;
          ++m_nActiveElements;
        } 
        else if (m_nActiveElements) {
          bool valid_set = m_elements[m_nElements].set(0,(*iter_boundaryLink),iter_PixelClusterColl,iter_PixelClusterCollEnd,theCluster,ctx);
          if(m_tools->isITkGeometry() && !valid_set) return false;
        } 
        else {
 
          continue;
        }
        m_elementsMap[m_nElements] = m_nElements;
        if (++m_nElements==300) break;
      }   
    }   
    else if (SCT) {
      InDet::SCT_ClusterCollection::const_iterator iter_stripClusterColl, iter_StripClusterCollEnd;
      const InDet::SCT_ClusterCollection *clustersOnElement = (*SCTc).indexFindPtr(id);

      if (clustersOnElement!=nullptr && clustersOnElement->begin()!=clustersOnElement->end()) {
 
        iter_stripClusterColl = clustersOnElement->begin();
        iter_StripClusterCollEnd = clustersOnElement->end();

        for (iter_cluster=lSiCluster.begin(); iter_cluster!=lSiCluster.end(); ++iter_cluster) {
          if ((*iter_cluster)->detectorElement()==detectorElement) {
            if (m_nclusters==0){
              m_firstElement = m_nElements;
            }
            else{
              m_lastElement  = m_nElements;
            }
            ++m_nclusters;
            m_ndfcut+=1;
            theCluster=(*iter_cluster);
            iter_cluster=lSiCluster.erase(iter_cluster);
            break;
          }
        }
        bool valid_set = m_elements[m_nElements].set(1,(*iter_boundaryLink),iter_stripClusterColl,iter_StripClusterCollEnd,theCluster,ctx);
        if(m_tools->isITkGeometry() && !valid_set) return false;
        ++m_nActiveElements;
      } 
      else if (m_nActiveElements) {
        bool valid_set = m_elements[m_nElements].set(0,(*iter_boundaryLink),iter_stripClusterColl,iter_StripClusterCollEnd,theCluster,ctx);
        if(m_tools->isITkGeometry() && !valid_set) return false;
      }
      else {
        continue;
      }
      m_elementsMap[m_nElements] = m_nElements;
      if (++m_nElements==300) break;
    } 

    if (m_firstElement == m_nElements-1) {
      up = m_nElements-1;         
      if (!m_elements[up].firstTrajectorElement(Tp, ctx)) return false;  
    }
    else if (theCluster) {

      if (!m_elements[m_nElements-1].ForwardPropagationWithoutSearch(m_elements[up],ctx)) {
        return false;
      }
      up = m_nElements-1;
      if (m_elements[m_nElements-1].xi2F() <= Xi2cut) {
        last=up;
      } 
      else {
        if (m_tools->isITkGeometry()) return false;
        else{
          ndfwrong+=m_elements[m_nElements-1].ndf();
          if (ndfwrong > 3) return false;
          m_ndfcut-=m_elements[m_nElements-1].ndf();
          --m_nclusters;
          m_elements[m_nElements-1].eraseClusterForwardPropagation();
        }
      }
    }
    else if (m_nclusters && !lSiCluster.empty()) {
      if (!m_elements[m_nElements-1].ForwardPropagationWithoutSearch(m_elements[up], ctx)) {
        if(not m_tools->isITkGeometry() and m_elements[m_nElements-1].cluster()) return false;
        --m_nElements;
        if (m_elements[m_nElements-1].detstatus()) --m_nActiveElements;
      } 
      else {
        if (m_elements[m_nElements-1].inside()<0) ++m_nholes;
        up = m_nElements-1;
      }
    } 
  } 

  if (!lSiCluster.empty()) {
    rquality = false;
    return false;
  }
  if (not m_tools->isITkGeometry() && ndfwrong && m_ndfcut < 6) return false;

  m_ndf = m_ndfcut;
  if (m_tools->isITkGeometry() || m_ndfcut > 6) m_ndfcut = 6;

  int n = m_nElements-1;
  for (; n>0; --n) {
    if (m_elements[n].detstatus()>=0) break;
  }
  m_nElements = n+1;
  
  for (; n>0; --n) {
    if (m_elements[n].detstatus() == 1) {
      m_elements[n].lastActive();
      break;
    }
  }

  int m         = m_firstElement+1;
  m_lastElement = last            ;
  for (n = m; n!=m_lastElement; ++n) {
    InDet::SiTrajectoryElement_xk& En = m_elements[m_elementsMap[n]];
    if (En.cluster() || En.inside() <= 0) m_elementsMap[m++] = m_elementsMap[n];
  }

  if (m!=n) {
    m_lastElement = m;
    for (; n!=m_nElements; ++n){
       m_elementsMap[m++] = m_elementsMap[n];
    }
    m_nElements = m;
  }
  if (!m_tools->bremNoise()) return true;

  for (n=m_lastElement; n!=m_nElements; ++n) {
    m_elements[m_elementsMap[n]].bremNoiseModel();
  }
  return true;
}

// Seacrh cluster compatible with track parameters
 
bool InDet::SiTrajectory_xk::trackParametersToClusters
(const InDet::PixelClusterContainer*                       PIXc      ,
 const InDet::SCT_ClusterContainer*                        SCTc      ,
 const Trk::TrackParameters                              & Tp        ,
 std::vector<const InDet::SiDetElementBoundaryLink_xk*>    & DE      ,
 std::multimap<const Trk::PrepRawData*,const Trk::Track*>& PT        ,
 std::vector<const InDet::SiCluster*>                      & lSiCluster, 
 const EventContext& ctx)
{
  m_nElements = 0;
  m_ndf       = 0;
 
  std::multimap<double,const InDet::SiCluster*> xi2cluster;
 
  std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator iter_boundaryLink,endBoundaryLinks=DE.end();
  std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator t, te =PT.end();
 
  double xi2Cut = .5;
  int    ndfCut =  6;
 
  for (iter_boundaryLink=DE.begin(); iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
 
    const InDetDD::SiDetectorElement* detectorElement = (*iter_boundaryLink)->detElement();
    IdentifierHash id = detectorElement->identifyHash();
 
    bool sct = detectorElement->isSCT();
 
    if (!sct) {
      InDet::PixelClusterCollection::const_iterator sib, sie;
      const InDet::PixelClusterCollection *w = (*PIXc).indexFindPtr(id);
 
      if (w!=nullptr && w->begin()!=w->end()) {
        sib = w->begin();
        sie = w->end  ();
      } else {
        continue;
      }
      if (!m_elements[0].ForwardPropagationForClusterSeach(m_nElements,Tp,(*iter_boundaryLink),sib,sie,ctx)) return false;
    } else {
      InDet::SCT_ClusterCollection::const_iterator sib, sie;
      const InDet::SCT_ClusterCollection *w = (*SCTc).indexFindPtr(id);
 
      if (w!=nullptr && w->begin()!=w->end()) {
        sib = w->begin();
        sie = w->end  ();
      } else {
        continue;
      }
      if (!m_elements[0].ForwardPropagationForClusterSeach(m_nElements,Tp,(*iter_boundaryLink),sib,sie,ctx)) return false;
    }
 
    for (int i=0; i!=m_elements[0].nlinksF(); ++i) {
    
      double x = m_elements[0].linkF(i).xi2();
 
      if (sct) {
        t = PT.find(m_elements[0].linkF(i).cluster());
        if (t!=te && (*t).second->measurementsOnTrack()->size() >= 10) continue;
      } else {
        x*=.5;
      }
 
      if (x <= xi2Cut) xi2cluster.insert(std::make_pair(x,m_elements[0].linkF(i).cluster()));
      break;
    }
    ++m_nElements;
  }
  
  if (xi2cluster.size() < 3) return false;
 
  std::multimap<double,const InDet::SiCluster*>::iterator xc = xi2cluster.begin(), xce = xi2cluster.end();
 
  for (; xc!=xce; ++xc) {
    lSiCluster.push_back((*xc).second);
    (*xc).second->detectorElement()->isSCT() ? m_ndf+=1 : m_ndf+=2;
    if ( m_ndf >= ndfCut ) break;
  }
 
  return m_ndf >= 6;
}

// Seacrh cluster compatible with global positions
 
bool InDet::SiTrajectory_xk::globalPositionsToClusters
(const InDet::PixelClusterContainer*                       PIXc      ,
 const InDet::SCT_ClusterContainer*                        SCTc      ,
 const std::vector<Amg::Vector3D>                          & Gp        ,
 std::vector<const InDet::SiDetElementBoundaryLink_xk*>    & DE        ,
 std::multimap<const Trk::PrepRawData*,const Trk::Track*>& PT        ,
 std::vector<const InDet::SiCluster*>                      & lSiCluster)
{
  std::vector<const InDet::SiDetElementBoundaryLink_xk*>::iterator iter_boundaryLink = DE.begin(), endBoundaryLinks = DE.end();
  std::vector<Amg::Vector3D>::const_iterator g,gb = Gp.begin(), ge = Gp.end();
  InDet::PixelClusterCollection::const_iterator pib, pie;
  InDet::SCT_ClusterCollection::const_iterator sib, sie;
  std::multimap<const Trk::PrepRawData*,const Trk::Track*>::const_iterator t, te =PT.end();
 
  Trk::PatternTrackParameters Tp;
  
  double pv[ 5]={0.,0.,0.,0.,0.};
  double cv[15]={ .1 ,
                  0. , .1,
                  0. , 0.,.001,
                  0. , 0.,  0.,.001,
                  0. , 0.,  0.,  0.,.00001};
 
  double xi2Cut = 10.;
  m_ndf         = 0  ;
 
  for (; iter_boundaryLink!=endBoundaryLinks; ++iter_boundaryLink) {
 
    const InDetDD::SiDetectorElement* d  = (*iter_boundaryLink)->detElement();
    IdentifierHash                    id = d->identifyHash ();
    const Trk::Surface*               su = &d->surface();
    const Trk::PlaneSurface* pla = static_cast<const Trk::PlaneSurface*>(su);
    if (!pla) continue;
 
    const Amg::Transform3D&  tr = pla->transform();
    double Ax[3] = {tr(0,0),tr(1,0),tr(2,0)};
    double Ay[3] = {tr(0,1),tr(1,1),tr(2,1)};
    double Az[3] = {tr(0,2),tr(1,2),tr(2,2)};
    double x0    = tr(0,3);
    double y0    = tr(1,3);
    double z0    = tr(2,3);
    double zcut  = .001  ;
    
    bool sct = d->isSCT();
    if (!sct) {
      const InDet::PixelClusterCollection *w = (*PIXc).indexFindPtr(id);
      if (w!=nullptr && w->begin()!=w->end()) {
        pib = w->begin();
        pie = w->end  ();
      } else {
        continue;
      }
    } else {
      zcut = 1.;
      const InDet::SCT_ClusterCollection *w = (*SCTc).indexFindPtr(id);
      if (w!=nullptr && w->begin()!=w->end()) {
        sib = w->begin();
        sie = w->end  ();
      } else {
        continue;
      }
    }
 
    for (g=gb; g!=ge; ++g) {
      
      double dx = (*g).x()-x0;
      double dy = (*g).y()-y0;
      double dz = (*g).z()-z0;
      double z  = dx*Az[0]+dy*Az[1]+dz*Az[2];
      if (std::abs(z) > zcut) continue;
 
      pv[0]     = dx*Ax[0]+dy*Ax[1]+dz*Ax[2];
      pv[1]     = dx*Ay[0]+dy*Ay[1]+dz*Ay[2];
      //setParametersWithCovariance detects whether su is 'owned' elsewhere
      //the ownership patterns could be improved, though.
      //coverity[MULTIPLE_INIT_SMART_PTRS]
      Tp.setParametersWithCovariance(su,pv,cv);
 
      if (!sct) m_elements[0].CloseClusterSeach(Tp, (*iter_boundaryLink), pib, pie);
      else      m_elements[0].CloseClusterSeach(Tp, (*iter_boundaryLink), sib, sie);
      const InDet::SiCluster* c =  m_elements[0].cluster();
      if (!c || m_elements[0].xi2F() > xi2Cut) continue;
      if (sct) {
        t = PT.find(c);
        if (t!=te && (*t).second->measurementsOnTrack()->size() >= 10) continue;
      }
      sct ? m_ndf+=1 : m_ndf+=2;
      lSiCluster.push_back(c);
    }
  }
  return m_ndf >= 6;
}

// Backward test initial trajectory
 
bool InDet::SiTrajectory_xk::backwardSmoother(bool TWO, const EventContext& ctx)
{
  if (m_firstElement >= m_lastElement) return false;
 
  // Trajectory difference test
  //
  int m = m_lastElement;
  for (; m>=m_firstElement; --m) {
    if (m_elements[m_elementsMap[m]].difference()) break;
  }
  if (m < m_firstElement) return true;
 
  if (!m_elements[m_elementsMap[m_lastElement]].lastTrajectorElement()) return false;
  
  int firstElement = m_lastElement       ;
  int maxholes     = m_tools->maxholes ();
  int maxdholes    = m_tools->maxdholes();
  m_nclustersNoAdd = 0                   ;
  m_difference     = 0                   ;
 
  m     = m_lastElement-1;
  int n = m              ;
 
  for (; m>=m_firstElement; --m) {
 
    InDet::SiTrajectoryElement_xk& En = m_elements[m_elementsMap[m+1]];
    InDet::SiTrajectoryElement_xk& Em = m_elements[m_elementsMap[m  ]];
 
    if (!Em.BackwardPropagationSmoother(En,TWO,ctx)) {
      
      if (m == m_firstElement) break;
 
      for (int i=m+1; i!=m_nElements; ++i) m_elementsMap[i-1] = m_elementsMap[i];
      --m_lastElement;
      --m_nElements;
      --firstElement;
      continue;
    }
 
    if ((Em.cluster() && Em.clusterOld()) && (Em.cluster()!=Em.clusterOld())) ++m_difference;
 
    if     (Em.cluster()) {
      firstElement  = m;
    }
    else                  {
      n=m;
      if (Em.clusterNoAdd()) ++m_nclustersNoAdd;
      if (Em.nholesB() > maxholes || Em.dholesB() > maxdholes) {
        ++m_difference; break;
      }
    }
  } 
 
  m_firstElement = firstElement                       ;
  m              = m_elementsMap[firstElement]        ;
  n              = m_elementsMap[     n      ]        ;
  m_nclusters    = m_elements[m].nclustersB()         ;
  m_nholes       = m_elements[m].nholesB   ()         ;
  m_nHolesBefore      = m_elements[n].nholesB   ()-m_nholes;
  m_ndf          = m_elements[m].ndfB()               ;
  if (m_ndf < m_ndfcut) return false;
 
  // Erase trajectory elements with big distance from the track
  //
  m = firstElement+1;
  n = m;
  for (; m!=m_lastElement; ++m) {
 
    InDet::SiTrajectoryElement_xk& Em = m_elements[m_elementsMap[m]];
    if (Em.inside() > 0) {
      if (Em.detstatus() > 0) --m_nActiveElements;
    }
    else                 {
      m_elementsMap[n++] = m_elementsMap[m];
    }
  }
  m_lastElement = n;
 
  // Test number trajector elemenst with clusters
  //
  if (m_nActiveElements < m_tools->clustersmin() && m_nholes+m_nHolesAfter) return false;
  if (n!=m) {
    for (; m!=m_nElements; ++m) m_elementsMap[n++] = m_elementsMap[m];
    m_nElements = n;
  }
  return true;
}

// Backward trajectory extension
 
bool InDet::SiTrajectory_xk::backwardExtension(int itmax, const EventContext& ctx)
{
  if (m_firstElement >= m_lastElement) return false;
  int L = m_firstElement;
  if (L==0) return true;
 
  int                     MPbest[300]         ;
  int                     TE    [100]         ;
  const InDet::SiCluster* CL    [100]         ;
  double                  XI2B  [100]         ;
  //Local variable PUB uses 27200 bytes of stack space
  //coverity[STACK_USE]
  Trk::PatternTrackParameters PUB[100]        ;
  Trk::PatternTrackParameters  PA             ;
 
  int    maxholes       = m_tools->maxholes ();
  int    maxdholes      = m_tools->maxdholes();
  const int itm         = itmax-1             ;
  int    it             = 0                   ;
  int    itbest         = 0                   ;
  int    qbest          =-100                 ;
  int    nbest          = 0                   ;
  int    ndfbest        = 0                   ;
  int    lbest          = L                   ;
  int    hbest          = 0                   ;
  int    hbestb         = 0                   ;
  int    nclbest        = 0                   ;
  int    ndcut          = 3                   ;
  int    F              = L                   ;
  double Xi2best        = 0.                  ;
 
  m_elements[m_elementsMap[F]].setNdist(0);
 
  for (; it!=itmax; ++it) {
 
    int  l = F;
    int  lastElementWithExpHit = F;
 
    for (--F; F>=0; --F) {
      
      InDet::SiTrajectoryElement_xk& El = m_elements[m_elementsMap[F+1]];
      InDet::SiTrajectoryElement_xk& Ef = m_elements[m_elementsMap[F  ]];
 
      if (!Ef.BackwardPropagationFilter(El, ctx))  break;
      
      if (Ef.cluster()) {
        lastElementWithExpHit = F;
        l = F;
      }
      else if (Ef.inside() < 0) {
        lastElementWithExpHit = F;
        if (Ef.nholesB() > maxholes || Ef.dholesB() > maxdholes) break;
      }
 
      int nm = Ef.nclustersB()+F;
      if (Ef.ndist() >  ndcut ||
      nm < nclbest ||
      (nm == nclbest && Ef.xi2totalB() > Xi2best) ) break;
    } 
 
    int    fl = F;
    if (fl<0) fl = 0;
 
    int    m  = m_elementsMap.at(l);
    int    nc = m_elements[m].nclustersB();
 
    if (it==0 && nc==m_nclusters) return true;
 
    int    np = m_elements[m].npixelsB();
    int    nh = m_elements[m].nholesB();
    int    nd = m_elements[m_elementsMap[fl]].ndist();
    int    q  = nc-nh;
    double X  = m_elements[m].xi2totalB();
 
    if     ( (q > qbest) || (q==qbest && X < Xi2best ) ) {
 
      qbest   = q                                          ;
      nbest   = 0                                          ;
      ndfbest = 0                                          ;
      hbest   = nh                                         ;
      hbestb  = m_elements[m_elementsMap[lastElementWithExpHit]].nholesB()-nh  ;
      itbest  = it                                         ;
      Xi2best = X                                          ;
      PA      = m_elements[m_elementsMap.at(l)].parametersUB();
 
      if (fl==0 && nd < ndcut) ndcut = nd;
 
      if (fl!=0 || nd > 0 || np < 3) {
 
        lbest = l-1;
        for (int i=l; i!=L; ++i) {
 
          InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap[i]];
   
          if (Ei.inside() <= 0 ) {
            MPbest[++lbest] = i;
            if (Ei.cluster()) {
          CL[nbest]   = Ei.cluster();
          XI2B[nbest] = Ei.xi2B();
          PUB[nbest]  = Ei.parametersUB();
          TE[nbest++] = lbest;
          ndfbest += Ei.ndf();
        }
          }
        }
      }
      else    {
 
        l     =   -1;
        lbest = fl-1;
        for (int i=fl; i!=L; ++i) {
 
          InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap[i]];
 
          if (Ei.inside() <= 0 && ++lbest >=0 ) {
            MPbest[lbest] = lbest;
            if (Ei.cluster()) {
          CL[nbest]   = Ei.cluster();
          XI2B[nbest] = Ei.xi2B();
          PUB[nbest]  = Ei.parametersUB();
          TE[nbest++] = lbest;
          ndfbest += Ei.ndf();
          if (l<0) l=lbest;
            }
            m_elementsMap.at(lbest) = m_elementsMap.at(i);
          }
   
        }
 
        int dn = L-1-lbest;
 
        if (dn!=0) {
 
          for (int i=L; i!= m_nElements; ++i) {
            m_elementsMap.at(i-dn)=m_elementsMap.at(i);
          }
 
          L            -=dn;
          m_nElements  -=dn;
          m_lastElement-=dn;
        }
      }
      nclbest = m_nclusters+nbest;
    }
    
    F = -1;
    if (l<=0) l=1;
    bool cl = false;
    double Xn = 0.;
 
    for (; l < L; ++l) {
      InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap.at(l)];
 
      if (Ei.cluster() && Ei.isNextClusterHoleB(cl,Xn))  {
        int nm = l+Ei.nclustersB();
        if (!cl) {
      if (Ei.dist() < -2. && Ei.ndist() > ndcut-1 ) continue;
      --nm;
    }
        if (nm < nclbest || (nm == nclbest && Xn > Xi2best)) continue;
        F=l; break;
      }
    }
 
    if (F < 0 ) break;
    if (it!=itm) if (!m_elements[m_elementsMap[F]].addNextClusterB()) break;
  }
  if (it == itmax) --it;
  if (!nbest) return true;
 
  m_nholes       = hbest ;
  m_nHolesBefore      = hbestb;
  m_nclusters   += nbest ;
  m_ndf         += ndfbest;
  m_firstElement = TE[0] ;
  m_elements[m_elementsMap[TE[0]]].setParametersB(PA);
 
  int dn = L-1-lbest;
 
  if (dn != 0) {
 
    m_nElements  -=dn;
    m_lastElement-=dn;
 
    int n = m_firstElement;
    for (; n <= lbest     ; ++n) m_elementsMap[n]=m_elementsMap[MPbest[n]];
    for (; n!= m_nElements; ++n) m_elementsMap[n]=m_elementsMap[n    +dn ];
  }
 
  if (itbest==it) return true;
  
  for (int n = L-1-dn; n>=0; --n) {
 
    InDet::SiTrajectoryElement_xk& En = m_elements[m_elementsMap[n]];
    int m = nbest-1;
    for (; m>=0; --m) if (TE[m]==n) break;
 
    if (m>=0) {
      if (m_tools->useFastTracking()) {
    En.setClusterB(CL[m],XI2B[m]);
    En.setParametersB(PUB[m]);
      }
      else En.setCluster(CL[m]);
      if (--nbest==0) break;
    }
    else     {
      if (m_tools->useFastTracking()) En.setClusterB(  nullptr  ,10000.);
      else En.setCluster(  nullptr  );
    }
  } 
  return true;
}

// Forward trajectory extension
 
bool InDet::SiTrajectory_xk::forwardExtension(bool smoother,int itmax, const EventContext& ctx)
{
  const double pi2 = 2.*M_PI;
  const double pi = M_PI;
  
  if (m_firstElement >= m_lastElement) return false;

  int extensionStartIndex = m_lastElement;
  int lastElementOnTraj = m_nElements-1;

  if (smoother) {
 
    extensionStartIndex = m_firstElement;

    if (m_elements[m_elementsMap[extensionStartIndex]].difference()) {

      if (!m_elements[m_elementsMap[extensionStartIndex]].firstTrajectorElement(false)) return false;
      for (++extensionStartIndex; extensionStartIndex<=m_lastElement; ++extensionStartIndex) {
        InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[extensionStartIndex-1]];
        InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[extensionStartIndex  ]];
        if (!thisElement.ForwardPropagationWithoutSearch(previousElement,ctx)) return false;
      }
    }
    else{
      bool diff = false;
      for (++extensionStartIndex; extensionStartIndex<=m_lastElement; ++extensionStartIndex) {
 
        InDet::SiTrajectoryElement_xk& previousElement = m_elements[m_elementsMap[extensionStartIndex-1]];
        InDet::SiTrajectoryElement_xk& thisElement = m_elements[m_elementsMap[extensionStartIndex  ]];
        if (!diff) {
          diff = thisElement.difference();
          if (diff) {
            if (!thisElement.addNextClusterF(previousElement,thisElement.cluster())) return false;
          }
        }
        else      {
          if (!thisElement.ForwardPropagationWithoutSearch(previousElement,ctx)) return false;
        }
      } 
    } 
    --extensionStartIndex;  
  } 

  if ( extensionStartIndex== lastElementOnTraj) return true;
 
  // Search best forward trajectory prolongation
  int                     MP    [300]                              ;
  int                     MPbest[300]                              ;
  int                     TE    [100]                              ;
  const InDet::SiCluster* CL    [100]                              ;

  int    maxholes       = m_tools->maxholes ()                     ;
  int    maxdholes      = m_tools->maxdholes()                     ;
  const int itm         = itmax-1                                  ;
  int    iteration      = 0                                        ;
  int    itbest         = 0                                        ;
  int    qbest          =-100                                      ;
  int    nbest          = 0                                        ;
  int    ndfbest        = 0                                        ;
  int    hbest          = 0                                        ;
  int    hbeste         = 0                                        ;
  int    ndbest         = 0                                        ;
  int    ndcut          = 3                                        ;
  int    nclbest        = 0                                        ;
  int    lbest          = extensionStartIndex                      ;
  int    index_currentElement              = extensionStartIndex   ;
  int    M              = extensionStartIndex                      ;
  MP    [M]             = extensionStartIndex                      ;
  double Xi2best        = 0.                                       ;
  const double dfmax    = 2.2                                      ;

  double f0             = m_elements[m_elementsMap[m_firstElement]].parametersUF().parameters()[2];
 
 
  m_elements[m_elementsMap.at(index_currentElement)].setNdist(0);

  for (; iteration!=itmax; ++iteration) {
 
    int  lastElementWithSeenHit = index_currentElement;
    int  lastElementWithExpHit  = index_currentElement;
    int  index_previousElement  = index_currentElement;
    int  mLastCluster = M;
    int  Cm = nclbest-lastElementOnTraj;
    bool haveHole = false;

    for (++index_currentElement; index_currentElement!=m_nElements; ++index_currentElement) {
      
      InDet::SiTrajectoryElement_xk& prevElement    = m_elements[m_elementsMap.at(index_previousElement)];
      InDet::SiTrajectoryElement_xk& currentElement = m_elements[m_elementsMap.at(index_currentElement )];

      if (!currentElement.ForwardPropagationWithSearch(prevElement, ctx)) {

        if (!currentElement.isBarrel() || index_previousElement!=index_currentElement-1) break;

        int index_auxElement = index_currentElement;
        for (; index_auxElement!=m_nElements; ++index_auxElement) {
          if (!m_elements[m_elementsMap[index_auxElement]].isBarrel() ) break;
        }
        if (index_auxElement==m_nElements) break;

        index_currentElement = index_auxElement-1;
        continue;
      } 

      else {    
        index_previousElement = index_currentElement;
      }  
      
      MP[++M] = index_currentElement;

      if (currentElement.cluster()) {
        if (not m_tools->isITkGeometry()) {
          double df = std::abs(currentElement.parametersUF().parameters()[2]-f0);
          if (df > pi) df = pi2-df;
          if (df > dfmax) break;
        }
        lastElementWithExpHit  = index_currentElement;
        lastElementWithSeenHit = index_currentElement;
        mLastCluster = M;
        haveHole = false;
      }

      else if (currentElement.inside() < 0  ) {
        lastElementWithExpHit=index_currentElement;
        if (currentElement.nholesF() > maxholes || currentElement.dholesF() > maxdholes) break;
 
        haveHole = true;
        if (not m_tools->isITkGeometry()) {
          double df = std::abs(currentElement.parametersPF().parameters()[2]-f0);
          if (df > pi) df = pi2-df;
          if (df > dfmax ) break;
        }
      }

      else if (not m_tools->isITkGeometry()) {
        double df = std::abs(currentElement.parametersPF().parameters()[2]-f0);
        if (df > pi) df = pi2-df;
        if (df > dfmax) break;
      }
      int nm = currentElement.nclustersF()-index_currentElement;
      if (     currentElement.ndist() > ndcut                       
           ||  nm < Cm                                              
           || (nm == Cm && currentElement.xi2totalF() > Xi2best)    
          ) break;
    } 


    int    m  = m_elementsMap[lastElementWithSeenHit];
    int    nc = m_elements[m].nclustersF();
    int    nh = m_elements[m].nholesF();
    m_nHolesAfter = m_elements[m_elementsMap[lastElementWithExpHit]].nholesF()-nh;

    if (iteration==0 && nc==m_nclusters) return true;

    int    lastElementProcessed = index_currentElement;
    if (lastElementProcessed==m_nElements) --lastElementProcessed;

    int    nd =  m_elements[m_elementsMap[lastElementProcessed]].ndist();
    int    q  = nc-nh;
    double X  = m_elements[m].xi2totalF();
    if     ( (q > qbest) || (q==qbest && X < Xi2best ) ) {
      qbest   = q;
      nbest   = 0;
      ndfbest = 0;
      hbest   = nh;
      hbeste  = m_elements[m_elementsMap[lastElementWithExpHit]].nholesF()-nh;
      itbest  = iteration;
      lbest   = extensionStartIndex ;
      Xi2best = X ;
      ndbest  = nd;
      if (lastElementProcessed==lastElementOnTraj && nd < ndcut) ndcut = nd;

      for (int j=extensionStartIndex+1; j<=mLastCluster; ++j) {
        int    i  = MP[j];
        InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap[i]];
        if (Ei.inside() <= 0) {
          MPbest[++lbest] = i;
          if (Ei.cluster()) {
            CL[nbest]   = Ei.cluster();
            TE[nbest++] = lbest;
            ndfbest += Ei.ndf();
          }
        }
      }
      nclbest = m_nclusters+nbest;
      if ( (nclbest >= 14 && !haveHole) || (lastElementProcessed==lastElementOnTraj && ndbest == 0)) break;
    } 

    index_currentElement = -1;
    bool cl = false;
    int nb = lastElementOnTraj-nclbest-1;
    double Xn;

    
    for (int j=mLastCluster; j!=extensionStartIndex; --j) {
      int i = MP[j];
      InDet::SiTrajectoryElement_xk& Ei = m_elements[m_elementsMap[i]];

      if (i!=lastElementOnTraj && Ei.cluster() && Ei.isNextClusterHoleF(cl,Xn)) {

        int nm = nb-i+Ei.nclustersF();

        if (!cl) {
          if (Ei.dist() < -2. && Ei.ndist() > ndcut-1) continue;
        }
        else  ++nm;
        
        if (nm < 0 || (nm == 0 &&   Xn > Xi2best)) continue;
        index_currentElement = i;
        M = j;
        break;
      }
    } 

    if (index_currentElement < 0 ) break;
    if (iteration!=itm && !m_elements.at(m_elementsMap.at(index_currentElement)).addNextClusterF()) break;
  }

  if (iteration == itmax) --iteration;

  if (!nbest) return true;

  m_nholes      = hbest      ;
  m_nHolesAfter     = hbeste     ;
  m_nclusters  += nbest      ;
  m_ndf        += ndfbest    ;
  m_lastElement = TE[nbest-1];
  m_nElements   = m_lastElement+1;

  if (m_lastElement != MPbest[m_lastElement]) {
    for (int n = extensionStartIndex+1; n<=m_lastElement; ++n){
       m_elementsMap[n]=m_elementsMap[MPbest[n]];
    }
  }
  if (itbest==iteration) return true;

  int mb =  0;
  index_currentElement      = -1;

  for (int n = extensionStartIndex+1; n!=m_nElements; ++n) {
    InDet::SiTrajectoryElement_xk& En = m_elements[m_elementsMap[n]];
 
    int m = mb;
    for (; m!=nbest; ++m) if (TE[m]==n) break;
    if (m!=nbest) {
      if (CL[m]!=En.cluster()) {
        if (index_currentElement<0) {
          index_currentElement=n;
          En.addNextClusterF(m_elements[m_elementsMap[n-1]],CL[m]);
        } 
        else    {     
          En.setCluster(CL[m]);
        }
      } 
      if (++mb == nbest) break;
    }  
    else {
      if (En.cluster()) {
        if (index_currentElement<0) {
          index_currentElement = n;
          En.addNextClusterF(m_elements[m_elementsMap[n-1]],nullptr);
        }
        else    {     
          En.setCluster(nullptr);
        }
      }
    }
  }  

  if (index_currentElement < 0 || m_lastElement == index_currentElement) {
    return true;
  }  

  for (++index_currentElement; index_currentElement<=m_lastElement; ++index_currentElement) {
    InDet::SiTrajectoryElement_xk& prevElement = m_elements[m_elementsMap.at(index_currentElement-1)];
    InDet::SiTrajectoryElement_xk& currentElement = m_elements[m_elementsMap.at(index_currentElement)];
    if (!currentElement.ForwardPropagationWithoutSearch(prevElement, ctx)) return false;
  }
  return true;
}

// Get clusters list from trajectory
 
void InDet::SiTrajectory_xk::getClusters
(std::vector<const InDet::SiCluster*>& Cl) const
{
  for (int i = m_firstElement; i<=m_lastElement; ++i) {
    int m = m_elementsMap[i];
    if (m_elements[m].cluster()) Cl.push_back(m_elements[m].cluster());
  }
}

// Forward filter without search
 
bool InDet::SiTrajectory_xk::forwardFilter(const EventContext& ctx)
{
  int L = m_firstElement;
  
  if (!m_elements[m_elementsMap[L]].firstTrajectorElement(false)) return false;
 
  for (++L; L<=m_lastElement; ++L) {
 
    InDet::SiTrajectoryElement_xk& El = m_elements[m_elementsMap[L-1]];
    InDet::SiTrajectoryElement_xk& Ef = m_elements[m_elementsMap[L  ]];
    if (!Ef.ForwardPropagationWithoutSearch(El,ctx)) return false;
  }
  return true;
}
 

// Filter with precise clusters error
 
bool InDet::SiTrajectory_xk::filterWithPreciseClustersError(const EventContext& ctx)
{
  int L       = m_firstElement;
  int I       = 0             ;
 
  if(!m_elements[m_elementsMap[L]].cluster()) return false;
  m_elementsMap[I] =  m_elementsMap[L];
 
  for(++L; L<=m_lastElement; ++L) {
 
    int K = m_elementsMap[L];
    if(m_elements[K].cluster() ||
       m_elements[K].clusterNoAdd() ||
       m_elements[K].inside() < 0 ) m_elementsMap[++I] = K;
  }
  m_firstElement = 0  ;
  m_lastElement  = I  ;
  m_nElements    = I+1;
 
  // Forward filter
  //
  L = 0;
 
  // firstTrajectorElement with correction = true
  if(!m_elements[m_elementsMap[L]].firstTrajectorElement(true)) return false;
 
  for(++L; L<=m_lastElement; ++L) {
 
    InDet::SiTrajectoryElement_xk& El = m_elements[m_elementsMap[L-1]];
    InDet::SiTrajectoryElement_xk& Ef = m_elements[m_elementsMap[L  ]];
 
    if(!Ef.ForwardPropagationWithoutSearchPreciseWithCorrection(El, ctx)) return false;
  }
 
  // Backward smoother
  //
  if(!m_elements[m_elementsMap[m_lastElement]].lastTrajectorElementPrecise()) return false;
 
  int m = m_lastElement-1;
  for(; m>=0; --m) {
 
    InDet::SiTrajectoryElement_xk& En = m_elements[m_elementsMap[m+1]];
    InDet::SiTrajectoryElement_xk& Em = m_elements[m_elementsMap[m  ]];
 
    if(!Em.BackwardPropagationPrecise(En, ctx)) return false;
  }
 
  return true;
}
 

// Test order of detector elements
 
bool InDet::SiTrajectory_xk::goodOrder()
{
  int    L    = m_firstElement  ;
  int    n    = m_elementsMap[L];
  double step = 0.              ;
  bool  order = true            ;
 
  for (++L; L<=m_lastElement; ++L) {
 
    int m = m_elementsMap[L];
    if (!m_elements[m].cluster() &&
    !m_elements[m].clusterNoAdd() &&
    m_elements[m].inside()>=0) continue;
 
    double stp =  m_elements[m].step(m_elements[n]);
    if     ( step   ==   0.) step = stp  ;
    else if ((step*stp) < 0.) {order = false; break;} 
    n = m;
  }
  if (order) return true;
 
  L = m_firstElement ;
  n = m_elementsMap[L];
  Amg::Vector3D gp = m_elements[n].globalPosition();
  double rad = gp.x()*gp.x()+gp.y()*gp.y();
  for (++L; L<=m_lastElement; ++L) {
  
    int m = m_elementsMap[L];
    if (!m_elements[m].cluster() &&
    !m_elements[m].clusterNoAdd() &&
    m_elements[m].inside()>=0) continue;
 
    gp = m_elements[m].globalPosition();
    double R = gp.x()*gp.x()+gp.y()*gp.y();
    if (R < rad) return false;
    rad = R;
    n=m;
  }
  return true;
}

// Sort of detector elements in step order
 
void  InDet::SiTrajectory_xk::sortStep()
{
  int    L  = m_firstElement;
  int    LA = m_firstElement;
 
  for (++L; L<=m_lastElement; ++L) {
    
    int m = m_elementsMap[L];
    if (!m_elements[m].cluster() &&
    !m_elements[m].clusterNoAdd() &&
    m_elements[m].inside()>=0) continue;
 
    m_elementsMap[++LA] = m;
  }
 
  m_lastElement = LA;
  L = m_firstElement;
  m_nElements = LA+1;
 
  bool   nc = true;
  bool   so = true;
  double ds = m_elements[m_elementsMap[LA]].step()-m_elements[m_elementsMap[L]].step();
  
  if (ds > 0.) {                // Sort in increase order
 
    while(nc) {
      nc = false;
      int m = L, n = L+1;
      for (; n<=LA; ++n) {
 
        int Mn = m_elementsMap[n];
        int Mm = m_elementsMap[m];
 
        if (m_elements[Mn].step() < m_elements[Mm].step()) {
          if (m_elements[Mn].step(m_elements[Mm]) < 0.) {
            m_elementsMap[m] = Mn;
            m_elementsMap[n] = Mm;
            nc = true; so = false;
          }
        }
        ++m;
      }
    }
  }
  else {
 
    while(nc) {                  // Sort in decrease order
      nc = false;
      int m = L, n = L+1;
      for (; n<=LA; ++n) {
 
        int Mn = m_elementsMap[n];
        int Mm = m_elementsMap[m];
 
        if (m_elements[Mn].step() > m_elements[Mm].step()) {
          if (m_elements[Mn].step(m_elements[Mm]) > 0.) {
            m_elementsMap[m] = Mn;
            m_elementsMap[n] = Mm;
            nc = true; so = false;
          }
        }
        ++m;
      }
    }
  }
  if (so) return;
 
  // Search first detector elements with cluster
  //
  int n = L;
  for (; n<= LA; ++n) {
    
    int e = m_elementsMap[n];
    if     (m_elements[e].cluster()) break;
    if     (m_elements[e].clusterNoAdd()) --m_nclustersNoAdd;
    else if (m_elements[e].inside() < 0 &&
         m_elements[e].detstatus()>=0) {--m_nholes; ++m_nHolesBefore;}
  }
 
  // Search last detector elements with cluster
  //
  int m = LA;
  for (; m>=n  ; --m) {
 
    int e = m_elementsMap[m];
    if     (m_elements[e].cluster()) break;
    if     (m_elements[e].clusterNoAdd()) --m_nclustersNoAdd;
    else if (m_elements[e].inside() < 0 &&
         m_elements[e].detstatus()>=0) {--m_nholes; ++m_nHolesAfter;}
  }
  m_firstElement = n;
  m_lastElement  = m;
 
}

// Test possibility for trajectory to jump through perigee 
 
bool InDet::SiTrajectory_xk::jumpThroughPerigee()
{
  int i = m_firstElement;
  double St = m_elements[m_elementsMap[m_lastElement]].step()-m_elements[m_elementsMap[i]].step();
 
  for (; i<=m_lastElement; ++i) {
    int m = m_elementsMap[i];
    if (m_elements[m].cluster() &&
    (m_elements[m].ndf()!=2 || (m_elements[m].stepToPerigee()*St) <= 0.)) break;
 
    if (m_elements[m].cluster()){
      --m_nclusters;
      m_ndf -= m_elements[m].ndf();
    }
    else if (m_elements[m].clusterNoAdd()) --m_nclustersNoAdd;
    else                                   --m_nholes        ;
  }
 
  if (i == m_firstElement) return false;
  m_firstElement = i;
  return true;
}

// Trajectory quality without optimization
 
double InDet::SiTrajectory_xk::quality() const
{
  int    holes   = 0 ;
  double quality = 0.;
 
  for (int i = m_firstElement; i<=m_lastElement; ++i) {
    quality+=m_elements[m_elementsMap[i]].quality(holes);
  }
  return quality;
}

// Trajectory quality with otimization
 
double InDet::SiTrajectory_xk::qualityOptimization() 
{
  int    lE = m_firstElement;
  int    h  = 0             ;
  double q  = 0             ;
  double qM = 0.            ;
 
  for (int i = m_firstElement; i<=m_lastElement; ++i) {
    int  m = m_elementsMap[i];
    q+=m_elements[m].quality(h);
    if (m_elements[m].cluster() && q > qM) {
      qM = q;
      lE = i;
    }   
  }
 
  if (lE == m_firstElement) return -100;
 
  int fE             = lE;
  int nclustersNoAdd = 0 ;
  int nclusters      = 0 ;
  int nholes         = 0 ;
  int dholes         = 0 ;
  int ndf            = 0 ;
  h                  = 0 ;
  q                  = 0.;
  qM                 = 0.;
  
  for (int i = lE; i>=m_firstElement; --i) {
 
    int  m = m_elementsMap[i];
    q+=m_elements[m].quality(h);
    
    if (m_elements[m].cluster()) {
 
      ++nclusters;
      ndf+=m_elements[m].ndf();
 
      if (q > qM) {
        qM = q;
    fE = i;
        m_nclusters      = nclusters     ;
        m_nclustersNoAdd = nclustersNoAdd;
        m_nholes         = nholes        ;
        m_dholes         = dholes        ;
        m_ndf            = ndf           ;
      }
 
    }
    else if (m_elements[m].clusterNoAdd()) {
      ++nclustersNoAdd;
    }
    else if (m_elements[m].inside() < 0 && m_elements[m].detstatus() >=0) {
      ++nholes;
      if (h > dholes) dholes = h;
    }
  }
 
  if (fE==lE || m_nclusters+m_nclustersNoAdd < m_tools->clustersmin()) return -100.;
  m_firstElement = fE;
  m_lastElement  = lE;
  return qM;
}

// Trajectory conversion to TrackStateOnSurface for next tracks
 
Trk::TrackStates
InDet::SiTrajectory_xk::convertToNextTrackStateOnSurface()
{
  int i=0;
  for (; i!=m_ntos; ++i) {
    if (m_itos[i]+1 < m_elements[m_atos[i]].ntsos()) {++m_itos[i]; break;}
    m_itos[i] = 0;
  }
  if (i==m_ntos) {
    return Trk::TrackStates();
  }
 
   auto dtsos = Trk::TrackStates();
 
  for (i=0; i!=m_ntos; ++i) {
 
    Trk::TrackStateOnSurface* tsos = m_elements[m_atos[i]].tsos(m_itos[i]);
    if (tsos) dtsos.push_back(tsos);
  }
  return dtsos;
}

// pT of the first elementtrajectory
 
double InDet::SiTrajectory_xk::pTfirst() const
{
  int n = m_firstElement  ; if (n <0 || n>=300) return 0.;
  n     = m_elementsMap[n]; if (n <0 || n>=300) return 0.;
  int s = m_elements[n].status();
  if (s<=1) return 0.;
  return m_elements[n].parametersUB().momentum().perp();
}

void InDet::SiTrajectory_xk::updateHoleSearchResult(){

  m_patternHoleOutcome = InDet::PatternHoleSearchOutcome{0,0,0,0,0,true};

  bool prevIsSctHole = false;

  
  for (int theEle  = m_firstElement+1; theEle<m_lastElement; ++theEle) {
    int m = m_elementsMap[theEle];
    InDet::SiTrajectoryElement_xk & theElement = m_elements[m];
    bool isPix = theElement.ndf() == 2;
    bool isSCTHole=false;

    if (theElement.cluster() || theElement.clusterNoAdd()) {
      prevIsSctHole = false;
      continue;
    } 
    
    else {
      std::unique_ptr<const Trk::TrackParameters> pars {theElement.trackParameters(true,0)};
      Trk::BoundaryCheckResult boundaryStatus = m_tools->boundaryCheckTool()->boundaryCheck(*pars);
      switch (boundaryStatus){
        case Trk::BoundaryCheckResult::Candidate: 
          if (isPix) {
            ++m_patternHoleOutcome.nPixelHoles;
          }
          else {
            ++m_patternHoleOutcome.nSCTHoles;
            isSCTHole=true;
          } 
          break;
        case Trk::BoundaryCheckResult::Insensitive: 
        case Trk::BoundaryCheckResult::OnEdge:      
        case Trk::BoundaryCheckResult::Outside:     
        case Trk::BoundaryCheckResult::Error:       
          break;
        case Trk::BoundaryCheckResult::DeadElement:
          if (isPix){
            ++m_patternHoleOutcome.nPixelDeads;
          }
          else{
            ++m_patternHoleOutcome.nSCTDeads;
          }
          break;
      } 
    } 

    if (isSCTHole && prevIsSctHole){
      prevIsSctHole = false;
      ++m_patternHoleOutcome.nSCTDoubleHoles;
    }
    else {
      prevIsSctHole = isSCTHole;
    }
  } 
  if (m_patternHoleOutcome.nPixelHoles+m_patternHoleOutcome.nSCTHoles > m_tools->maxholes()
      || m_patternHoleOutcome.nSCTDoubleHoles > m_tools->maxdholes()){
     m_patternHoleOutcome.passPatternHoleCut = false;
  }
}