Skip to content
Snippets Groups Projects

Fixed an error in the X projection fit

Merged Fixed an error in the X projection fit
lohenry-fixXFit into master
Merged Louis Henry requested to merge lohenry-fixXFit into master
Compare
and
1 file
+ 28
− 42
Compare changes
  • Side-by-side
  • Inline
Pr/PrAlgorithms/src/PrHybridSeeding.cpp
+ 28
− 42
@@ -606,7 +606,7 @@ bool PrHybridSeeding::addStereo<CASE>::fitTrackSimultaneouslyXY( Pr::Hybrid::See
const float dxdy = hit.dxDy();
const float yOnTrack = track.yOnTrack( hit );
const float dz = 0.001f * ( hit.z( yOnTrack ) - Pr::Hybrid::zReference );
const float deta = dz * dz * ( 1.f + dz * track.dRatio() );
const float deta = dz * dz * ( 0.001f + dz * track.dRatio() );
const float wdz = w * dz;
const float weta = w * deta;
const float wdxdy = w * dxdy;
@@ -638,7 +638,7 @@ bool PrHybridSeeding::addStereo<CASE>::fitTrackSimultaneouslyXY( Pr::Hybrid::See
if ( !decomp ) return false;
decomp.Solve( rhs );
rhs[1] *= 1.e-3f;
rhs[2] *= 1.e-6f;
rhs[2] *= 1.e-9f;
rhs[4] *= 1.e-3f;
rhs[3] -= rhs[4] * Pr::Hybrid::zReference;
track.updateParameters( rhs[0], rhs[1], rhs[2], rhs[3], rhs[4] );
@@ -796,55 +796,41 @@ bool PrHybridSeeding::fitXProjection( Pr::Hybrid::SeedTrackX& track, unsigned in
//---LoH: called 30k times
// if ( track.size() < m_minXPlanes ) return false;
float mat[6];
std::fill( mat, mat + 3, 0.f );
std::fill( mat, mat + 6, 0.f );
float rhs[3];
// First loop to set some things
std::fill( rhs, rhs + 3, 0.f );
const float dRatio = track.dRatio();
for ( const auto& modHit : track.hits() ) {
const PrHit& hit = *( modHit.hit );
float w = hit.w(); // squared
const float dz = 0.001f * ( hit.z() - Pr::Hybrid::zReference );
const PrHit& hit = *( modHit.hit );
float w = hit.w(); // squared
const float dz = 0.001f * ( hit.z() - Pr::Hybrid::zReference );
const float deta = dz * dz * ( 0.001f + dRatio * dz );
const float dist = track.distance( hit );
mat[0] += w;
mat[1] += w * dz;
mat[2] += w * dz * dz;
mat[3] += w * deta;
mat[4] += w * dz * deta;
mat[5] += w * deta * deta;
rhs[0] += w * dist;
rhs[1] += w * dist * dz;
rhs[2] += w * dist * deta;
}
for ( unsigned int loop = 0; 3 > loop; ++loop ) {
std::fill( mat + 3, mat + 6, 0.f );
std::fill( rhs, rhs + 3, 0.f );
for ( const auto& modHit : track.hits() ) {
// for( auto itH = track.hits().begin(); track.hits().end() != itH; ++itH ){
const PrHit& hit = *( modHit.hit );
const float dRatio = track.dRatio();
float w = hit.w(); // squared
const float dz = 0.001f * ( hit.z() - Pr::Hybrid::zReference );
const float deta = dz * dz * ( 1.f + dRatio * dz );
const float dist = track.distance( hit );
mat[3] += w * deta;
mat[4] += w * dz * deta;
mat[5] += w * deta * deta;
rhs[0] += w * dist;
rhs[1] += w * dist * dz;
rhs[2] += w * dist * deta;
}
ROOT::Math::CholeskyDecomp<float, 3> decomp( mat ); //---LoH: can probably be made more rapidly
if ( !decomp ) { return false; }
// Solve linear system
decomp.Solve( rhs );
rhs[1] *= 1.e-3f;
rhs[2] *= 1.e-6f;
// protect against unreasonable track parameter corrections
//---LoH: commented as it slows down more that accelerates things
// if ( std::fabs( rhs[0] ) > 1.e4f || std::fabs( rhs[1] ) > 5.f || std::fabs( rhs[2] ) > 1.e-3f ) return false;
// Small corrections
track.updateParameters( rhs[0], rhs[1], rhs[2] );
// Put back later faster maybe
if ( loop > 0 && std::abs( rhs[0] ) < 5e-5f && std::abs( rhs[1] ) < 5e-8f && std::abs( rhs[2] ) < 5e-11f ) {
break;
}
}
ROOT::Math::CholeskyDecomp<float, 3> decomp( mat ); //---LoH: can probably be made more rapidly
if ( !decomp ) { return false; }
// Solve linear system
decomp.Solve( rhs );
rhs[1] *= 1.e-3f;
rhs[2] *= 1.e-9f;
// protect against unreasonable track parameter corrections
//---LoH: commented as it slows down more that accelerates things
// if ( std::fabs( rhs[0] ) > 1.e4f || std::fabs( rhs[1] ) > 5.f || std::fabs( rhs[2] ) > 1.e-3f ) return false;
// Small corrections
track.updateParameters( rhs[0], rhs[1], rhs[2] );
// Compute some values on the track
auto hitChi2 = [&track]( const ModPrHit& hit ) { return track.chi2( *hit.hit ); };
const auto chi2s = track.hits() | ranges::views::transform( hitChi2 );
float sum( 0 );
float sum( 0.f );
for ( auto chi2 : chi2s ) sum += chi2;
track.setChi2( sum, track.hits().size() - 3 );
return ranges::max( chi2s ) < m_maxChi2HitsX[iCase];