diff --git a/Phys/VertexFit/src/ParticleVertexFitter.cpp b/Phys/VertexFit/src/ParticleVertexFitter.cpp index e350181cce45c7f9825cd120ba454f75b19fbec3..6bddbab7a10255076237515b233e9d7d7fe8ae3e 100644 --- a/Phys/VertexFit/src/ParticleVertexFitter.cpp +++ b/Phys/VertexFit/src/ParticleVertexFitter.cpp @@ -78,12 +78,12 @@ private: StatusCode fit( const LHCb::Particle::ConstVector&, LHCb::Vertex&, LHCb::Particle*, const IGeometryInfo& ) const; private: - Gaudi::Property<int> m_maxnumiter{this, "MaxNumIter", 5}; - Gaudi::Property<double> m_maxdchisq{this, "MaxDeltaChi2", 0.01}; - Gaudi::Property<bool> m_extrapolateTracksWithoutVelo{this, "extrapolateTracksWithoutVelo", false}; - PublicToolHandle<ITrackStateProvider> m_stateprovider{this, "StateProvider", "TrackStateProvider"}; - ToolHandle<ITrajPoca> m_trajpoca{"TrajPoca"}; - const LHCb::IParticlePropertySvc* m_ppsvc{nullptr}; + Gaudi::Property<int> m_maxnumiter{this, "MaxNumIter", 5}; + Gaudi::Property<double> m_maxdchisq{this, "MaxDeltaChi2", 0.01}; + Gaudi::Property<bool> m_extrapolateTtracks{this, "extrapolateTtracks", false}; + PublicToolHandle<ITrackStateProvider> m_stateprovider{this, "StateProvider", "TrackStateProvider"}; + ToolHandle<ITrajPoca> m_trajpoca{"TrajPoca"}; + const LHCb::IParticlePropertySvc* m_ppsvc{nullptr}; mutable Gaudi::Accumulators::MsgCounter<MSG::WARNING> m_matrix_inversion_failed{this, "Problem inverting matrix!", 0}; }; @@ -591,7 +591,8 @@ StatusCode ParticleVertexFitter::fit( const LHCb::Particle::ConstVector& origdau // 1. if there are 'resonance' daughters, take the vertex position of the first. // 2. if not, use the states of the first two tracks with velo or composites // 3. if not, try with downstream tracks and trajpoca - bool posinitialized{false}; + bool posinitialized{false}; + Gaudi::XYZPoint position; if ( counttypes[Resonance] > 0 ) { // try 1 @@ -630,7 +631,45 @@ StatusCode ParticleVertexFitter::fit( const LHCb::Particle::ConstVector& origdau trajs[ntrajs++] = &( ownedtrajs.back() ); } } - if ( ntrajs == 2 ) { + bool parallelTracks = false; + if ( ntrajs == 2 && m_extrapolateTtracks ) { + // this is NOT default behaviour + // For T-tracks it does not make sense to use 3-D POCA starting position + // Trajectory in YZ plane is approx straight line + // calculate the yz intersection z and extrapolate there + + const float ty0 = daughters[0]->proto()->track()->firstState().ty(); + const float ty1 = daughters[1]->proto()->track()->firstState().ty(); + + if ( ty0 == ty1 ) { + debug() << "tracks are parallel in yz plane. Cannot calculate yz intersection. Default to 3-D POCA starting " + "position. " + << endmsg; + parallelTracks = true; + } else { + const float z0 = daughters[0]->proto()->track()->firstState().z(); + const float z1 = daughters[1]->proto()->track()->firstState().z(); + + const float y0 = daughters[0]->proto()->track()->firstState().y(); + const float y1 = daughters[1]->proto()->track()->firstState().y(); + + // calculate the y-axis intersection of each track + const float c0 = y0 - ty0 * z0; + const float c1 = y1 - ty1 * z1; + + // calculate the point in the yz-plane where the tracks intersect + const float yzIntersectionZ = ( c1 - c0 ) / ( ty0 - ty1 ); + const float yzIntersectionY = ty0 * yzIntersectionZ + c1; + + position.SetZ( yzIntersectionZ ); + position.SetY( yzIntersectionY ); + position.SetX( 0 ); + posinitialized = true; + } + } + + if ( ( ntrajs == 2 && !m_extrapolateTtracks ) || ( ntrajs == 2 && parallelTracks ) ) { + // this is the default behaviour double mu0( 0 ), mu1( 0 ); Gaudi::XYZVector deltaX; StatusCode sc = m_trajpoca->minimize( *( trajs[0] ), mu0, *( trajs[1] ), mu1, deltaX, 0.1 * Gaudi::Units::mm ); @@ -675,7 +714,7 @@ StatusCode ParticleVertexFitter::fit( const LHCb::Particle::ConstVector& origdau } } break; case TrackWithoutVelo: { - if ( !m_extrapolateTracksWithoutVelo ) { + if ( !m_extrapolateTtracks ) { // this is the default behaviour const LHCb::Track* track = p->proto()->track(); const LHCb::TrackTraj* tracktraj = m_stateprovider->trajectory( *( p->proto()->track() ), geometry );