Trk::FitProcedure Class Reference

#include <FitProcedure.h>

Collaboration diagram for Trk::FitProcedure:

Classes
struct	Cache

Public Member Functions
	FitProcedure (bool constrainedAlignmentEffects, bool extendedDebug, bool lineFit, ToolHandle< IIntersector > &rungeKuttaIntersector, ToolHandle< IIntersector > &solenoidalIntersector, ToolHandle< IIntersector > &straightLineIntersector, const ToolHandle< IPropagator > &stepPropagator, const Volume *indetVolume=0, int maxIterations=25, int useStepPropagator=0)
	~FitProcedure (void)=default
const FitProcedureQuality &	execute (FitProcedure::Cache &cache, bool asymmetricCaloEnergy, MsgStream &log, std::vector< FitMeasurement * > &measurements, FitParameters *&parameters, const FitQuality *perigeeQuality=0, bool for_iPatTrack=false) const
void	setMinIterations (int minIter)
bool	constrainedAlignmentEffects () const

Static Public Member Functions
static void	clear (FitProcedure::Cache &cache)
static Track *	constructTrack (FitProcedure::Cache &cache, const std::vector< FitMeasurement * > &measurements, FitParameters &parameters, const TrackInfo &trackInfo, const DataVector< const TrackStateOnSurface > *leadingTSOS=nullptr)
static Amg::MatrixX *	fullCovariance ()

Private Member Functions
	FitProcedure (const FitProcedure &)=delete
FitProcedure &	operator= (const FitProcedure &)=delete
	FitProcedure (FitProcedure &&)=delete
FitProcedure &	operator= (FitProcedure &&)=delete
void	calculateChiSq (FitProcedure::Cache &cache, std::vector< FitMeasurement * > &measurements) const
const ToolHandle< IIntersector > &	chooseIntersector (std::vector< FitMeasurement * > &measurements, const FitParameters &parameters) const

Static Private Member Functions
static void	reportQuality (FitProcedure::Cache &cache, const std::vector< FitMeasurement * > &measurements, const FitParameters &parameters)

Private Attributes
bool	m_constrainedAlignmentEffects
bool	m_extendedDebug
double	m_extremeOneOverP
const Trk::Volume *	m_indetVolume
double	m_largeRadius
bool	m_lineFit
int	m_maxIter
int	m_minIter
double	m_minPt
ToolHandle< IIntersector > &	m_rungeKuttaIntersector
const ToolHandle< IIntersector > &	m_solenoidalIntersector
ToolHandle< IIntersector > &	m_straightLineIntersector
const ToolHandle< IPropagator > &	m_stepPropagator
int	m_useStepPropagator

Detailed Description

Definition at line 44 of file FitProcedure.h.

Constructor & Destructor Documentation

FitProcedure() [1/3]

Trk::FitProcedure::FitProcedure	(	bool	constrainedAlignmentEffects,
		bool	extendedDebug,
		bool	lineFit,
		ToolHandle< IIntersector > &	rungeKuttaIntersector,
		ToolHandle< IIntersector > &	solenoidalIntersector,
		ToolHandle< IIntersector > &	straightLineIntersector,
		const ToolHandle< IPropagator > &	stepPropagator,
		const Volume *	indetVolume = 0,
		int	maxIterations = 25,
		int	useStepPropagator = 0 )

Definition at line 48 of file FitProcedure.cxx.

    : m_constrainedAlignmentEffects(constrainedAlignmentEffects),
      m_extendedDebug(extendedDebug),
      m_extremeOneOverP(1. / (10. * Gaudi::Units::TeV)),
      m_indetVolume(indetVolume),
      m_largeRadius(1000. * Gaudi::Units::mm),
      m_lineFit(lineFit),
      m_maxIter(maxIterations),
      m_minIter(0),
      m_minPt(0.05 * Gaudi::Units::GeV),
      m_rungeKuttaIntersector(rungeKuttaIntersector),
      m_solenoidalIntersector(solenoidalIntersector),
      m_straightLineIntersector(straightLineIntersector),
      m_stepPropagator(stepPropagator),
      m_useStepPropagator(useStepPropagator) {}

~FitProcedure()

Trk::FitProcedure::~FitProcedure ( void )

default

FitProcedure() [2/3]

Trk::FitProcedure::FitProcedure ( const FitProcedure & )

privatedelete

FitProcedure() [3/3]

Trk::FitProcedure::FitProcedure ( FitProcedure && )

privatedelete

Member Function Documentation

calculateChiSq()

void Trk::FitProcedure::calculateChiSq ( FitProcedure::Cache & cache, std::vector< FitMeasurement * > & measurements ) const

private

Definition at line 658 of file FitProcedure.cxx.

                                                    {
  // convergence criterion
  const double dChisqConv = 0.025;
 
  // compute total chisquared and sum of hit differences
  // flag hit with highest chisquared contribution (on entry if RoadFit)
  cache.chiSq = 0.;
  double driftResidual = 0.;
  double DSqMax = 0.;
  for (auto* m : measurements) {
    if (!m->numberDoF())
      continue;
    // if (m->isPerigee())
    // {
    //     cache.chiSq += cache.fitMatrices->perigeeChiSquared();
    //     continue;
    // }
 
    double residual = m->residual();
    double DiffSq = residual * residual;
    cache.chiSq += DiffSq;
    if (m->isPositionMeasurement()) {
      if (m->isDrift())
        driftResidual += residual;
      if (DiffSq > DSqMax) {
        DSqMax = DiffSq;
        cache.worstMeasurement = m->hitIndex() + 1;
        cache.chiSqWorst = std::min(999., DSqMax);
      }
    }
    if (m->is2Dimensional()) {
      residual = m->residual2();
      DiffSq = residual * residual;
      cache.chiSq += DiffSq;
      if (m->isPositionMeasurement()) {
        if (DiffSq > DSqMax) {
          DSqMax = DiffSq;
          cache.worstMeasurement = m->hitIndex() + 1;
          cache.chiSqWorst = std::min(999., DSqMax);
        }
      }
    }
  }
 
  // assess chi squared per degree of freedom (and its stability)
  if (cache.fitMatrices->numberDoF() > 0)
    cache.chiSq /= static_cast<double>(cache.fitMatrices->numberDoF());
  if (cache.iteration == 0) {
    cache.cutTaken = 0;
    cache.chRatio1 = 0.;
    cache.chRatio2 = 0.;
    cache.chiSqMin = cache.chiSq;
  }
 
  cache.chiSqOld = cache.chiSqMin;
  const double DChiSq = cache.chiSqOld - cache.chiSq;
  if (DChiSq > -dChisqConv) {
    cache.chiSqMin = cache.chiSq;
    cache.nCuts = 0;
  }
  if (cache.iteration > 0) {
    cache.chRatio2 = cache.chRatio1;
    cache.chRatio1 = cache.chiSq / cache.chiSqOld;
  }
  if (cache.fitMatrices->numberDriftCircles()) {
    cache.driftSumLast = cache.driftSum;
    cache.driftSum =
        driftResidual /
        static_cast<double>(cache.fitMatrices->numberDriftCircles());
  }
 
  //
  // debugging info
  if (cache.verbose) {
    *cache.log << "----------------------------------" << std::endl
               << std::setiosflags(std::ios::fixed)
               << " Debugging Info in ChiSquare method" << std::endl
               << " # of track-fit iterations " << std::setw(3)
               << cache.iteration << std::endl
               << " fit ChiSquared (per degree of freedom) " << std::setw(13)
               << std::setprecision(3) << cache.chiSq
               << "   # of degrees of freedom "
               << cache.fitMatrices->numberDoF() << std::endl
               << " ChSq Ratio1/2 " << std::setw(9) << std::setprecision(3)
               << cache.chRatio1 << std::setw(10) << std::setprecision(3)
               << cache.chRatio2 << std::endl
               << " driftResidual " << std::setw(9) << std::setprecision(3)
               << cache.driftSum << " #driftCircles "
               << cache.fitMatrices->numberDriftCircles() << std::endl
               << " CutTaken " << cache.cutTaken << std::endl
               << "----------------------------------" << std::endl
               << "   ";
 
    Trk::FitMeasurement::printHeading(*cache.log);
    int n = 0;
    for (auto* m : measurements) {
      *cache.log << std::setiosflags(std::ios::fixed) << std::setw(3) << ++n;
      if (m->isPerigee()) {
        *cache.log << " perigee     ";
        *cache.log << std::endl;
      } else {
        m->print(*cache.log);
      }
    }
  }
 
  //
  // check for possible convergence (nearConvergence forces extra iteration)
  if (!cache.cutStep && !cache.nCuts &&
      (cache.chiSq < 0.1 ||
       (cache.chRatio2 < 1.1 &&
        (std::abs(DChiSq) < dChisqConv ||
         std::abs((cache.chiSq - cache.chiSqOld) / cache.chiSqOld) < 0.01)))) {
    if ((cache.chiSq < 2.0 || cache.nearConvergence || cache.iteration == 1) &&
        cache.iteration >= m_minIter) {
      cache.convergence = true;
    } else {
      cache.nearConvergence = true;
      if (cache.verbose)
        *cache.log << MSG::VERBOSE << " near convergence " << endmsg;
    }
  } else {
    cache.nearConvergence = false;
  }
 
  // else take cutstep if divergent or oscillating
  cache.cutStep = false;
}

chooseIntersector()

const ToolHandle< IIntersector > & Trk::FitProcedure::chooseIntersector ( std::vector< FitMeasurement * > & measurements, const FitParameters & parameters ) const

private

Definition at line 789 of file FitProcedure.cxx.

                                           {
  if (m_lineFit) {
    return m_straightLineIntersector;
  }
 
  // decide which intersector to use for curved tracks (default RungeKutta)
  // ToolHandle<IIntersector>& intersector = m_rungeKuttaIntersector;
 
  // solenoidal intersector must start close to origin with last measurement
  // inside valid region
  for (std::vector<FitMeasurement*>::reverse_iterator m = measurements.rbegin();
       m != measurements.rend(); ++m) {
    if (!(**m).isPositionMeasurement())
      continue;
    if (!m_solenoidalIntersector->isValid(parameters.position(),
                                          (**m).position()))
      break;
    return m_solenoidalIntersector;
  }
 
  return m_rungeKuttaIntersector;
}

clear()

void Trk::FitProcedure::clear ( FitProcedure::Cache & cache )

static

Definition at line 73 of file FitProcedure.cxx.

                                                 {
  cache.fitMatrices->releaseMemory();
}

constrainedAlignmentEffects()

bool Trk::FitProcedure::constrainedAlignmentEffects ( ) const

inline

Definition at line 143 of file FitProcedure.h.

                                           {
    return m_constrainedAlignmentEffects;
  }

constructTrack()

Track * Trk::FitProcedure::constructTrack ( FitProcedure::Cache & cache, const std::vector< FitMeasurement * > & measurements, FitParameters & parameters, const TrackInfo & trackInfo, const DataVector< const TrackStateOnSurface > * leadingTSOS = nullptr )

static

Definition at line 77 of file FitProcedure.cxx.

                                       {
  // debug
  if (cache.debug) {
    reportQuality(cache, measurements, parameters);
  }
 
  // NB keep first and last measurements distinct i.e. separate TSOS (no
  // scatterers etc) NB trackParameters outwards from TSOS i.e. always last
  // FitMeas on surface
 
  // create vector of TSOS - reserve upper limit for size (+1 as starts with
  // perigee)
  auto trackStateOnSurfaces = std::make_unique<Trk::TrackStates>();
  unsigned size = measurements.size() + 1;
  if (leadingTSOS)
    size += leadingTSOS->size();
  trackStateOnSurfaces->reserve(size);
  std::unique_ptr<AlignmentEffectsOnTrack> alignmentEffects{};
  const FitMeasurement* fitMeasurement = measurements.front();
  const FitQualityOnSurface fitQoS{};
  std::unique_ptr<MaterialEffectsBase> materialEffects{};
  std::unique_ptr<MeasurementBase> measurementBase{};
  const Surface* surface = nullptr;
  std::unique_ptr<TrackParameters> trackParameters{};
  std::bitset<TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
      const defaultPattern;
  std::bitset<TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
      typePattern = defaultPattern;
 
  // start with (measured) perigee
  unsigned scatter = 0;
  unsigned tsos = 0;
  std::unique_ptr<Perigee> perigee(parameters.perigee());
  typePattern.set(TrackStateOnSurface::Perigee);
  trackStateOnSurfaces->push_back(new TrackStateOnSurface(
      fitQoS, std::move(measurementBase), std::move(perigee),
      std::move(materialEffects), typePattern, std::move(alignmentEffects)));
  ++tsos;
 
  // append leading TSOS to perigee
  if (leadingTSOS) {
    for (const auto* t : *leadingTSOS) {
      if (!(*t).type(Trk::TrackStateOnSurface::Perigee)) {
        trackStateOnSurfaces->push_back((*t).clone());
        ++tsos;
      }
    }
  }
 
  // then append the fitted TSOS
  for (auto* m : measurements) {
    if (m->isMaterialDelimiter())
      continue;
 
    // push back previous TSOS when fresh surface reached
    if (m->surface() != surface || alignmentEffects || m->alignmentEffects()) {
      if (surface) {
        if (typePattern == defaultPattern) {
          if (cache.debug) {
            *cache.log << MSG::DEBUG << " skip empty TSOS# " << tsos + 1;
            if (m->materialEffects())
              *cache.log << " with material";
            m->print(*cache.log);
            *cache.log << endmsg;
          }
        } else {
          // get the MeasuredParameters (with covariance)
          bool const withCovariance = true;
          trackParameters.reset(parameters.trackParameters(
              *cache.log, *fitMeasurement, withCovariance));
 
          if (!trackParameters) {
            *cache.log
                << MSG::WARNING
                << " fail track with incomplete return TSOS: no trackParameters"
                << endmsg;
            return nullptr;
          }
          typePattern.set(TrackStateOnSurface::Parameter);
          trackStateOnSurfaces->push_back(new TrackStateOnSurface(
              fitQoS, std::move(measurementBase), std::move(trackParameters),
              std::move(materialEffects), typePattern,
              std::move(alignmentEffects)));
          ++tsos;
        }
      }
      fitMeasurement = m;
      surface = m->surface();
      measurementBase.reset();
      materialEffects.reset();
      typePattern = defaultPattern;
      alignmentEffects.reset();
    } else {
      fitMeasurement = m;
      if (cache.verbose)
        *cache.log << MSG::VERBOSE << " tsos# " << tsos << " shared surface"
                   << endmsg;
    }
 
    // it's a measurement
    if (m->measurementBase()) {
      // create an extra TSOS if there is already a measurement on this surface
      // (dirty fix for pseudoMeasurements)
      if (measurementBase) {
        // get the MeasuredParameters (with covariance)
        bool const withCovariance = true;
        trackParameters.reset(parameters.trackParameters(
            *cache.log, *fitMeasurement, withCovariance));
        if (!trackParameters) {
          *cache.log
              << MSG::WARNING
              << " fail track with incomplete return TSOS: no trackParameters"
              << endmsg;
          return nullptr;
        }
        typePattern.set(TrackStateOnSurface::Parameter);
        trackStateOnSurfaces->push_back(new TrackStateOnSurface(
            fitQoS, std::move(measurementBase), std::move(trackParameters),
            std::move(materialEffects), typePattern,
            std::move(alignmentEffects)));
        ++tsos;
        fitMeasurement = m;
        materialEffects.reset();
        typePattern = defaultPattern;
        alignmentEffects.reset();
      }
 
      measurementBase = m->measurementBase()->uniqueClone();
      typePattern.set(TrackStateOnSurface::Measurement);
      if (m->isOutlier())
        typePattern.set(TrackStateOnSurface::Outlier);
    }
 
    // it's a CaloDeposit or Scatterer (scatterers may be fitted or not fitted)
    if (m->materialEffects()) {
      // update momentum to account for energy loss
 
      if (m->isEnergyDeposit()) {
        materialEffects = m->materialEffects()->uniqueClone();
        typePattern.set(TrackStateOnSurface::CaloDeposit);
      } else if (m->isScatterer()) {
        // set materialPattern as the scattering parameters are fitted
        std::bitset<MaterialEffectsBase::NumberOfMaterialEffectsTypes>
            typeMaterial;
        typeMaterial.set(MaterialEffectsBase::FittedMaterialEffects);
        const MaterialEffectsOnTrack* meot =
            dynamic_cast<const MaterialEffectsOnTrack*>(m->materialEffects());
        if (meot && meot->energyLoss())  // standard scatterer
        {
          auto energyLoss =
              std::unique_ptr<EnergyLoss>(meot->energyLoss()->clone());
          typeMaterial.set(Trk::MaterialEffectsBase::EnergyLossEffects);
          if (m->numberDoF())  // fitted scatterer
          {
            materialEffects = std::make_unique<MaterialEffectsOnTrack>(
                m->materialEffects()->thicknessInX0(),
                parameters.scatteringAngles(*m, scatter), std::move(energyLoss),
                m->materialEffects()->associatedSurface(), typeMaterial);
            ++scatter;
          } else  // unfitted (leading material)
          {
            materialEffects = std::make_unique<MaterialEffectsOnTrack>(
                m->materialEffects()->thicknessInX0(),
                parameters.scatteringAngles(*m), std::move(energyLoss),
                m->materialEffects()->associatedSurface(), typeMaterial);
          }
        } else  // no meot for special calo scattering centres
        {
          if (m->numberDoF())  // fitted scatterer
          {
            materialEffects = std::make_unique<MaterialEffectsOnTrack>(
                m->materialEffects()->thicknessInX0(),
                parameters.scatteringAngles(*m, scatter),
                m->materialEffects()->associatedSurface(), typeMaterial);
            ++scatter;
          } else  // unfitted (leading material)
          {
            materialEffects = std::make_unique<MaterialEffectsOnTrack>(
                m->materialEffects()->thicknessInX0(),
                parameters.scatteringAngles(*m),
                m->materialEffects()->associatedSurface(), typeMaterial);
          }
        }
 
        typePattern.set(TrackStateOnSurface::Scatterer);
      } else {
        *cache.log << MSG::WARNING
                   << " deprecated TrackStateOnSurface::InertMaterial"
                   << endmsg;
        materialEffects = m->materialEffects()->uniqueClone();
        typePattern.set(TrackStateOnSurface::InertMaterial);
      }
    }
 
    // additional perigee (e.g. at MS entrance)
    if (m->isPerigee()) {
      typePattern.set(TrackStateOnSurface::Perigee);
    }
 
    // or alignment effects
    else if (m->alignmentEffects()) {
      const AlignmentEffectsOnTrack& AEOT = *m->alignmentEffects();
      unsigned const align = m->alignmentParameter() - 1;
 
      *cache.log << MSG::VERBOSE << " Fitprocedure AEOT input deltaTranslation "
                 << AEOT.deltaTranslation() << " deltaAngle "
                 << AEOT.deltaAngle() << " output Trans "
                 << parameters.alignmentOffset(align) << " deltaAngle "
                 << parameters.alignmentAngle(align) << endmsg;
      alignmentEffects = std::make_unique<Trk::AlignmentEffectsOnTrack>(
          parameters.alignmentOffset(align), AEOT.sigmaDeltaTranslation(),
          parameters.alignmentAngle(align), AEOT.sigmaDeltaAngle(),
          AEOT.vectorOfAffectedTSOS(), *m->surface());
      typePattern.set(TrackStateOnSurface::Alignment);
    }
 
    // passive types: hole for now
    else if (m->isPassive()) {
      if (m->type() == hole)
        typePattern.set(TrackStateOnSurface::Hole);
    }
  }
 
  // remember the final TSOS !
  bool const withCovariance = true;
  trackParameters.reset(
      parameters.trackParameters(*cache.log, *fitMeasurement, withCovariance));
  if (!trackParameters) {
    *cache.log << MSG::WARNING
               << " fail track with incomplete return TSOS: no trackParameters"
               << endmsg;
    return nullptr;
  }
  typePattern.set(TrackStateOnSurface::Parameter);
  trackStateOnSurfaces->push_back(new TrackStateOnSurface(
      fitQoS, std::move(measurementBase), std::move(trackParameters),
      std::move(materialEffects), typePattern, std::move(alignmentEffects)));
  ++tsos;
 
  // construct track
  const double chiSquared = cache.chiSq * static_cast<double>(cache.numberDoF);
  Track* track = new Track(trackInfo, std::move(trackStateOnSurfaces),
                           std::make_unique<FitQuality>(chiSquared, cache.numberDoF));
 
  if (cache.verbose)
    *cache.log << MSG::VERBOSE << " track with " << tsos << " TSOS " << endmsg;
  return track;
}

execute()

const FitProcedureQuality & Trk::FitProcedure::execute	(	FitProcedure::Cache &	cache,
		bool	asymmetricCaloEnergy,
		MsgStream &	log,
		std::vector< FitMeasurement * > &	measurements,
		FitParameters *&	parameters,
		const FitQuality *	perigeeQuality = 0,
		bool	for_iPatTrack = false ) const

Definition at line 330 of file FitProcedure.cxx.

                                                                {
  // report start value
  cache.log = &log;
  if (cache.log->level() > MSG::DEBUG) {
    cache.debug = false;
    cache.verbose = false;
  } else {
    cache.debug = true;
    if (cache.log->level() < MSG::DEBUG)
      cache.verbose = true;
    *cache.log << MSG::DEBUG << "parameter start values:  ";
    parameters->print(*cache.log);
    *cache.log << endmsg;
  }
 
  // choose appropriate intersector
  const ToolHandle<IIntersector>& intersector =
      chooseIntersector(measurements, *parameters);
 
  // resize matrices
  int fitCode = cache.fitMatrices->setDimensions(measurements, parameters);
  if (fitCode) {
    cache.fitQuality = std::make_unique<FitProcedureQuality>(
        fitCode, cache.fitMatrices->numberDoF());
    if (cache.debug)
      reportQuality(cache, measurements, *parameters);
    return *cache.fitQuality;
  }
 
  // remaining initialization
  cache.chiSq = 0.;
  cache.chiSqWorst = 0.;
  cache.driftSum = 0.;
  cache.driftSumLast = 0.;
  cache.fitProbability = 0.;
  cache.iteration = -1;
  cache.numberDoF = cache.fitMatrices->numberDoF();
  cache.numberParameters = parameters->numberParameters();
  cache.worstMeasurement = 0;
  MeasurementProcessor measurementProcessor(
      asymmetricCaloEnergy, cache.fitMatrices->derivativeMatrix(), intersector,
      measurements, parameters, m_rungeKuttaIntersector, m_stepPropagator,
      m_useStepPropagator);
 
  // perigee or vertex used as measurements in fit
  if (measurements.front()->isPerigee()) {
    cache.fitMatrices->usePerigee(*measurements.front());
  }
 
  // set requested options and initial values
  const double ptInvCut = 1. / m_minPt;  // protection against trapped particles
  cache.cutStep = true;
  cache.convergence = false;
  cache.nearConvergence = false;
 
  // keep best (original if not reasonable quality) results
  double bestChiSq = cache.chiSqCut;
  std::optional<FitParameters> bestParameters = std::nullopt;
 
  // iteration loop to fit track parameters
  while (!fitCode && !cache.convergence) {
    bool forceIteration = false;
    if (cache.iteration > m_maxIter && bestParameters && !for_iPatTrack) {
      parameters->reset(*bestParameters);
      cache.convergence = true;
      cache.cutStep = false;
      if (cache.verbose)
        *cache.log << MSG::VERBOSE
                   << " convergence problem: accept after max iter " << endmsg;
    } else if (!cache.cutStep) {
      //  solve equations and update parameters
      if (!cache.iteration) {
        cache.fitMatrices->refinePointers();
        if (m_extendedDebug)
          cache.fitMatrices->checkPointers(*cache.log);
        if (cache.verbose)
          cache.fitMatrices->printDerivativeMatrix();
      }
 
      if (!cache.fitMatrices->solveEquations()) {
        fitCode = 11;  // fails matrix inversion
      } else if (parameters->fitEnergyDeposit() &&
                 !parameters->extremeMomentum() &&
                 std::abs(parameters->qOverP()) < m_extremeOneOverP) {
        if (cache.debug)
          *cache.log << MSG::DEBUG << " extremeMomentum " << endmsg;
        parameters->extremeMomentum(true);
        fitCode = cache.fitMatrices->setDimensions(measurements, parameters);
        bestChiSq = cache.chiSqCut;
        bestParameters = std::nullopt;
        forceIteration = true;
        cache.chiSq = 0.;
        cache.chiSqWorst = 0.;
        cache.driftSum = 0.;
        cache.driftSumLast = 0.;
        cache.numberParameters = parameters->numberParameters();
      }
      if (cache.verbose && !cache.iteration)
        cache.fitMatrices->printWeightMatrix();
    }
    ++cache.iteration;
 
    // report parameters
    if (cache.verbose) {
      *cache.log << MSG::VERBOSE << " ===== start iteration "
                 << cache.iteration;
      if (cache.iteration) {
        if (cache.cutStep)
          *cache.log << " ====== cutStep";
      } else {
        if (for_iPatTrack)
          *cache.log << "  ====== for_iPatTrack ";
      }
      parameters->printVerbose(*cache.log);
    }
 
    //  check for some error conditions (if none found yet)
    if (fitCode) {
      // e.g. fitCode == 11 (singular matrix)
    } else if (std::abs(parameters->ptInv0()) > ptInvCut) {
      fitCode = 8;  //  fail with pt below cutoff
    } else if (measurements.front()->isVertex() && m_indetVolume &&
               !m_indetVolume->inside(parameters->position())) {
      fitCode = 9;  // fail if vertex outside indet
    } else if (cache.iteration &&
               (std::abs(parameters->difference(3)) > 1.0 ||
                std::abs(parameters->difference(0)) > m_largeRadius) &&
               !measurements.front()->isVertex()) {
      if (std::abs(parameters->difference(3)) > 1.0) {
        fitCode = 10;  // fail with ill-defined cot_theta
      } else {
        fitCode = 9;  // fail crazy impact parameter
      }
    } else if (!fitCode) {
      //  extrapolate to each measurement and calculate derivatives
      if (!measurementProcessor.calculateFittedTrajectory(cache.iteration) ||
          !measurementProcessor.calculateDerivatives()) {
        fitCode = 5;  //  fail as trapped
        cache.fitQuality = std::make_unique<FitProcedureQuality>(
            cache.chiSq, cache.chiSqWorst, cache.fitProbability, fitCode,
            cache.iteration, parameters->numberAlignments(),
            cache.fitMatrices->numberDoF(), parameters->numberScatterers(),
            cache.worstMeasurement);
 
        if (cache.debug) {
          if (cache.verbose)
            *cache.log << endmsg;
          reportQuality(cache, measurements, *parameters);
        }
        return *cache.fitQuality;
      }
 
      //  have extrapolation and derivatives, calculate residual
      measurementProcessor.calculateResiduals();
 
      // check for remaining error conditions. If OK then compute chisquared.
      if (cache.iteration > m_maxIter && !cache.cutStep && for_iPatTrack) {
        fitCode = 6;  //  fail with no convergence
      } else if (cache.iteration == 4 && cache.chiSq > 1000. && for_iPatTrack) {
        fitCode = 7;  //  fail with too high chisquared
      } else if (!fitCode) {
        calculateChiSq(cache, measurements);
 
        // check for cutstep conditions if no significant chi2 improvement
        if (!forceIteration && !cache.convergence && cache.chRatio1 > 0.9) {
          double cutStep = 0.;
          if (cache.iteration > 4 &&
              cache.driftSum * cache.driftSumLast < -1.) {
            cache.cutStep = true;
            cutStep = std::abs(cache.driftSumLast /
                               (cache.driftSum - cache.driftSumLast));
            if (cutStep < 0.001)
              cutStep = 0.001;
            if (cache.verbose)
              *cache.log
                  << MSG::VERBOSE
                  << " take cutStep following chi2 increase on iteration "
                  << cache.iteration << "  chi2Ratio " << cache.chRatio1
                  << "   driftSum " << cache.driftSum << " prev "
                  << cache.driftSumLast << "   " << cutStep << endmsg;
          } else if (parameters->numberOscillations() > 2) {
            cache.cutStep = true;
            cutStep = 0.5;
            if (cache.verbose)
              *cache.log << MSG::VERBOSE
                         << " take cutStep as oscillating, iteration "
                         << cache.iteration << ", numberOscillations "
                         << parameters->numberOscillations() << endmsg;
          }
 
          // perform cutstep
          if (cache.cutStep) {
            cache.convergence = false;
            parameters->performCutStep(cutStep);
            if (cache.verbose)
              parameters->printVerbose(*cache.log);
            if (measurementProcessor.calculateFittedTrajectory(
                    cache.iteration)) {
              // note: derivatives should not be recalculated for cutstep
              measurementProcessor.calculateResiduals();
              calculateChiSq(cache, measurements);
              if (cache.verbose)
                *cache.log << "   after cutStep: "
                           << "  chi2Ratio " << cache.chRatio1 << "   driftSum "
                           << cache.driftSum << endmsg;
            }
          }
        }
 
        // keep current best parameters
        if (!bestParameters || cache.chiSq < bestChiSq) {
          bestChiSq = cache.chiSq;
          bestParameters = *parameters;
          parameters->resetOscillations();
        }
 
        if (bestParameters &&
            ((cache.convergence && cache.chiSq > bestChiSq + 0.5) ||
             (parameters->phiInstability() && cache.iteration == m_maxIter))) {
          parameters->reset(*bestParameters);
          if (cache.verbose) {
            *cache.log << MSG::VERBOSE << " revert to bestParameters ";
            parameters->printVerbose(*cache.log);
          }
          if (measurementProcessor.calculateFittedTrajectory(cache.iteration)) {
            measurementProcessor.calculateDerivatives();
            measurementProcessor.calculateResiduals();
            calculateChiSq(cache, measurements);
            cache.convergence = true;
          }
        }
 
        if (forceIteration)
          cache.convergence = false;
      }
    }  // if (std::abs(ptInv0) > ptInvCut)
    if (cache.verbose)
      *cache.log << endmsg;
 
    // try to rescue phi instability failures
    if (fitCode && cache.iteration && bestParameters &&
        !parameters->phiInstability() &&
        (**(measurements.rbegin())).position().perp() > m_largeRadius) {
      if (cache.verbose)
        *cache.log << MSG::VERBOSE << " phi instability " << endmsg;
      parameters->reset(*bestParameters);
      parameters->setPhiInstability();
      cache.cutStep = true;
      cache.convergence = false;
      fitCode = 0;
      cache.iteration = 0;
    }
  }  // while
 
  // store successful fit :
  if (!fitCode) {
    // set covariance in parameters class after inclusion of uncertainty from
    // field integral
    const Amg::MatrixX* fullCovariance = cache.fitMatrices->fullCovariance();
    if (fullCovariance) {
      Amg::MatrixX& finalCovariance = cache.fitMatrices->finalCovariance();
      if (!for_iPatTrack) {
        if (!m_lineFit)
          measurementProcessor.fieldIntegralUncertainty(*cache.log,
                                                        finalCovariance);
        measurementProcessor.propagationDerivatives();
      }
      parameters->covariance(&finalCovariance, fullCovariance);
 
      // fit quality
      if (perigeeQuality) {
        // take care when mixing normalized with unnormalized
        cache.chiSq = perigeeQuality->chiSquared() +
                      cache.chiSq * static_cast<double>(cache.numberDoF);
        cache.numberDoF += perigeeQuality->numberDoF();
        cache.chiSq /= static_cast<double>(cache.numberDoF);
      }
 
      // probability of chisquared
      cache.fitProbability = 1.;
      if (cache.numberDoF > 0 && cache.chiSq > 0.) {
        if (cache.chiSq < 100.) {
          const double chiSquared =
              cache.chiSq * static_cast<double>(cache.numberDoF);
          cache.fitProbability -=
              Genfun::CumulativeChiSquare(cache.numberDoF)(chiSquared);
        } else {
          cache.fitProbability = 0.;
        }
      }
 
      if (cache.verbose) {
        *cache.log << MSG::VERBOSE << " fit converged";
        if (cache.chiSqWorst > 6.25)
          *cache.log << " with possible outlier #" << cache.worstMeasurement
                     << " (residual " << std::sqrt(cache.chiSqWorst) << ")";
        *cache.log << endmsg;
      }
    } else {
      fitCode = 11;  //  singular weight matrix
    }
  }
 
  cache.fitQuality = std::make_unique<FitProcedureQuality>(
      cache.chiSq, cache.chiSqWorst, cache.fitProbability, fitCode,
      cache.iteration, parameters->numberAlignments(), cache.numberDoF,
      parameters->numberScatterers(), cache.worstMeasurement);
  if (cache.debug && (for_iPatTrack || fitCode))
    reportQuality(cache, measurements, *parameters);
 
  return *cache.fitQuality;
}

fullCovariance()

Amg::MatrixX * Trk::FitProcedure::fullCovariance ( )

static

Definition at line 646 of file FitProcedure.cxx.

                                         {
  // note const_cast - ughhh
  // return const_cast<Amg::MatrixX*>(cache.fitMatrices->fullCovariance());
  return nullptr;  // NOT mig5
}

operator=() [1/2]

FitProcedure & Trk::FitProcedure::operator= ( const FitProcedure & )

privatedelete

operator=() [2/2]

FitProcedure & Trk::FitProcedure::operator= ( FitProcedure && )

privatedelete

reportQuality()

void Trk::FitProcedure::reportQuality ( FitProcedure::Cache & cache, const std::vector< FitMeasurement * > & measurements, const FitParameters & parameters )

staticprivate

Definition at line 814 of file FitProcedure.cxx.

                                     {
  if (!cache.fitQuality)
    return;
 
  const int fitCode = cache.fitQuality->fitCode();
  if (fitCode) {
    *cache.log << MSG::DEBUG << "failure: fitCode " << fitCode;
    const std::string msg = "";
    switch (fitCode) {
      case 1:
        *cache.log << "  missing Trk::Surface ";
        break;
      case 2:
        *cache.log << "  too many measurements for fit matrix size: "
                   << measurements.size();
        break;
      case 3:
        *cache.log << "  too many parameters for fit matrix size: "
                   << parameters.numberParameters();
        break;
      case 4:
        *cache.log << "  unconstrained fit: negative numberDoF "
                   << cache.fitQuality->numberDoF();
        break;
      case 5:
        *cache.log << "  trapped in magnetic field: pT = "
                   << 1. / (parameters.ptInv0() * Gaudi::Units::GeV)
                   << "  at iteration# " << cache.fitQuality->iterations();
        break;
      case 6:
        *cache.log << "  no convergence: chiSq = " << cache.fitQuality->chiSq()
                   << "  at iteration# " << cache.fitQuality->iterations();
        break;
      case 7:
        *cache.log << "  enormous chi squared: chiSq = "
                   << cache.fitQuality->chiSq() << "  at iteration# "
                   << cache.fitQuality->iterations();
        break;
      case 8:
        *cache.log << "  below pT cutoff. pT = "
                   << 1. / (parameters.ptInv0() * Gaudi::Units::GeV)
                   << "  at iteration# " << cache.fitQuality->iterations();
        break;
      case 9:
        *cache.log << "  ill-defined impact parameter " << parameters.d0()
                   << "  with difference " << parameters.difference(0)
                   << "  at iteration# " << cache.fitQuality->iterations();
        break;
      case 10:
        *cache.log << "  ill-defined cotTheta " << parameters.cotTheta()
                   << "  with difference " << parameters.difference(3)
                   << "  at iteration# " << cache.fitQuality->iterations();
        break;
      case 11:
        *cache.log << "  singular matrix fails inversion:"
                   << "  at iteration# " << cache.fitQuality->iterations();
        break;
      case 12:
        *cache.log << "  maximum of one calorimeter permitted";
        break;
      case 13:
        *cache.log << "  NO derivativeMatrix available";
        break;
      default:
        break;
    };
    *cache.log << std::endl << endmsg;
  } else {
    *cache.log << MSG::DEBUG << "fitted parameter values: ";
    parameters.print(*cache.log);
    *cache.log << endmsg;
    *cache.log << MSG::DEBUG;
    cache.fitQuality->print(*cache.log);
    parameters.printCovariance(*cache.log);
    *cache.log << endmsg;
  }
}

setMinIterations()

void Trk::FitProcedure::setMinIterations

(

int

minIter

)

Definition at line 652 of file FitProcedure.cxx.

                                               {
  m_minIter = minIter;
  if (m_minIter > m_maxIter)
    m_maxIter = m_minIter;
}

Member Data Documentation

m_constrainedAlignmentEffects

bool Trk::FitProcedure::m_constrainedAlignmentEffects

private

Definition at line 165 of file FitProcedure.h.

m_extendedDebug

bool Trk::FitProcedure::m_extendedDebug

private

Definition at line 166 of file FitProcedure.h.

m_extremeOneOverP

double Trk::FitProcedure::m_extremeOneOverP

private

Definition at line 167 of file FitProcedure.h.

m_indetVolume

const Trk::Volume* Trk::FitProcedure::m_indetVolume

private

Definition at line 168 of file FitProcedure.h.

m_largeRadius

double Trk::FitProcedure::m_largeRadius

private

Definition at line 169 of file FitProcedure.h.

m_lineFit

bool Trk::FitProcedure::m_lineFit

private

Definition at line 170 of file FitProcedure.h.

m_maxIter

int Trk::FitProcedure::m_maxIter

private

Definition at line 171 of file FitProcedure.h.

m_minIter

int Trk::FitProcedure::m_minIter

private

Definition at line 172 of file FitProcedure.h.

m_minPt

double Trk::FitProcedure::m_minPt

private

Definition at line 173 of file FitProcedure.h.

m_rungeKuttaIntersector

ToolHandle<IIntersector>& Trk::FitProcedure::m_rungeKuttaIntersector

private

Definition at line 174 of file FitProcedure.h.

m_solenoidalIntersector

const ToolHandle<IIntersector>& Trk::FitProcedure::m_solenoidalIntersector

private

Definition at line 175 of file FitProcedure.h.

m_stepPropagator

const ToolHandle<IPropagator>& Trk::FitProcedure::m_stepPropagator

private

Definition at line 177 of file FitProcedure.h.

m_straightLineIntersector

ToolHandle<IIntersector>& Trk::FitProcedure::m_straightLineIntersector

private

Definition at line 176 of file FitProcedure.h.

m_useStepPropagator

int Trk::FitProcedure::m_useStepPropagator

private

Definition at line 178 of file FitProcedure.h.

---

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

• FitProcedure.h
• FitProcedure.cxx