|
|
|
#include <DualKalmanTrajectory.h>
Public Member Functions | |
| virtual DualKalmanTrajectory * | clone () const |
| DualKalmanTrajectory (const Trajectory::DataContainer &trajMeasurements, const TrajectoryStateOnSurface &referenceTsos, const std::vector< unsigned int > &forwardRecHitNums, const std::vector< unsigned int > &backwardRecHitNums, const MagneticField *magField, MaterialEffects materialEffects, PropagationDirection propDir, double mass, bool useBeamSpot, const reco::BeamSpot &beamSpot, int residualMethod) | |
| virtual | ~DualKalmanTrajectory () |
Protected Member Functions | |
| virtual bool | construct (const Trajectory::DataContainer &trajMeasurements, const TrajectoryStateOnSurface &referenceTsos, const std::vector< unsigned int > &forwardRecHitNums, const std::vector< unsigned int > &backwardRecHitNums, double mass, MaterialEffects materialEffects, const PropagationDirection propDir, const MagneticField *magField, bool useBeamSpot, const reco::BeamSpot &beamSpot, int residualMethod) |
| calculate members | |
| virtual ReferenceTrajectory * | construct (const Trajectory::DataContainer &trajMeasurements, const TrajectoryStateOnSurface &referenceTsos, const std::vector< unsigned int > &recHits, double mass, MaterialEffects materialEffects, const PropagationDirection propDir, const MagneticField *magField, bool useBeamSpot, const reco::BeamSpot &beamSpot) const |
| Method to get a single ReferenceTrajectory for a half of the trajectory. | |
| DualKalmanTrajectory (unsigned int nPar=0, unsigned int nHits=0) | |
| virtual AlgebraicVector | extractParameters (const TrajectoryStateOnSurface &referenceTsos) const |
| bool | fillKalmanPart (const Trajectory::DataContainer &trajMeasurements, const std::vector< unsigned int > &recHitNums, bool startFirst, unsigned int iNextHit, int residualMethod) |
| Fill that part of data members that is different from DualReferenceTrajectory. | |
| TrajectoryStateOnSurface | fillMeasurementAndError1 (const TransientTrackingRecHit::ConstRecHitPointer &hitPtr, unsigned int iHit, const TrajectoryMeasurement &trajMeasurement) |
| helper for 'unbiased residual' method (i.e. returns merged fwd/bwd states) | |
| TrajectoryStateOnSurface | fillMeasurementAndError2 (const TransientTrackingRecHit::ConstRecHitPointer &hitPtr, unsigned int iHit, const TrajectoryMeasurement &trajMeasurement, bool doPull) |
| void | fillTrajectoryPositions (const AlgebraicMatrix &projection, const TrajectoryStateOnSurface &tsos, unsigned int iHit) |
| fill trajectoryPositions | |
| LocalError | hitErrorWithAPE (const TransientTrackingRecHit::ConstRecHitPointer &hitPtr) const |
| local error including APE if APE is on | |
| const PropagationDirection | oppositeDirection (const PropagationDirection propDir) const |
A trajectory that provides derivatives w.r.t. to the track parameters as the DualReferenceTrajectory, i.e. extrapolating a helix with five parameters from the hit in the middle. But measurements, trajectory positions and uncertainties are taken from the Kalman track fit. More precisely, the uncertainties are those of the residuals, i.e. of measurements() - trajectoryPositions().
Currently three methods to set residual and error can be choosen via cfg: 0: the hitError approach, i.e. residual w.r.t. updated state, but error from hit 1: the unbiased residal approach as in validation 2: the pull approach, i.e. residal w.r.t. to updated state withits error
(last update by
)
(last update by
)
Definition at line 34 of file DualKalmanTrajectory.h.
| DualKalmanTrajectory::DualKalmanTrajectory | ( | const Trajectory::DataContainer & | trajMeasurements, |
| const TrajectoryStateOnSurface & | referenceTsos, | ||
| const std::vector< unsigned int > & | forwardRecHitNums, | ||
| const std::vector< unsigned int > & | backwardRecHitNums, | ||
| const MagneticField * | magField, | ||
| MaterialEffects | materialEffects, | ||
| PropagationDirection | propDir, | ||
| double | mass, | ||
| bool | useBeamSpot, | ||
| const reco::BeamSpot & | beamSpot, | ||
| int | residualMethod | ||
| ) |
Definition at line 33 of file DualKalmanTrajectory.cc.
References construct(), and ReferenceTrajectoryBase::theValidityFlag.
Referenced by clone().
: ReferenceTrajectoryBase(referenceTsos.localParameters().mixedFormatVector().kSize, forwardRecHitNums.size() + backwardRecHitNums.size() - 1, 0, 0) { theValidityFlag = this->construct(trajMeasurements, referenceTsos, forwardRecHitNums, backwardRecHitNums, mass, materialEffects, propDir, magField, useBeamSpot, beamSpot, residualMethod); }
| virtual DualKalmanTrajectory::~DualKalmanTrajectory | ( | ) | [inline, virtual] |
Definition at line 51 of file DualKalmanTrajectory.h.
{}
| DualKalmanTrajectory::DualKalmanTrajectory | ( | unsigned int | nPar = 0, |
| unsigned int | nHits = 0 |
||
| ) | [protected] |
Definition at line 55 of file DualKalmanTrajectory.cc.
: ReferenceTrajectoryBase(nPar, nHits, 0, 0) {}
| virtual DualKalmanTrajectory* DualKalmanTrajectory::clone | ( | void | ) | const [inline, virtual] |
Implements ReferenceTrajectoryBase.
Definition at line 53 of file DualKalmanTrajectory.h.
References DualKalmanTrajectory().
{ return new DualKalmanTrajectory(*this); }
| bool DualKalmanTrajectory::construct | ( | const Trajectory::DataContainer & | trajMeasurements, |
| const TrajectoryStateOnSurface & | referenceTsos, | ||
| const std::vector< unsigned int > & | forwardRecHitNums, | ||
| const std::vector< unsigned int > & | backwardRecHitNums, | ||
| double | mass, | ||
| MaterialEffects | materialEffects, | ||
| const PropagationDirection | propDir, | ||
| const MagneticField * | magField, | ||
| bool | useBeamSpot, | ||
| const reco::BeamSpot & | beamSpot, | ||
| int | residualMethod | ||
| ) | [protected, virtual] |
calculate members
Definition at line 61 of file DualKalmanTrajectory.cc.
References extractParameters(), fillKalmanPart(), TrajectoryStateOnSurface::hasError(), TrajectoryStateOnSurface::localError(), LocalTrajectoryError::matrix(), ReferenceTrajectoryBase::nMeasPerHit, oppositeDirection(), ReferenceTrajectoryBase::theDerivatives, ReferenceTrajectoryBase::theParameters, ReferenceTrajectoryBase::theRecHits, and ReferenceTrajectoryBase::theTrajectoryPositionCov.
Referenced by DualKalmanTrajectory().
{
ReferenceTrajectoryPtr fwdTraj = this->construct(trajMeasurements, refTsos, forwardRecHitNums,
mass, materialEffects,
propDir, magField, useBeamSpot, beamSpot);
ReferenceTrajectoryPtr bwdTraj = this->construct(trajMeasurements, refTsos, backwardRecHitNums,
mass, materialEffects,
this->oppositeDirection(propDir), magField,
useBeamSpot, beamSpot);
// edm::LogError("Alignment") << "@SUB=DualKalmanTrajectory::construct"
// << "valid fwd/bwd " << fwdTraj->isValid()
// << "/" << bwdTraj->isValid();
if (!fwdTraj->isValid() || !bwdTraj->isValid()) {
return false;
}
//
// Combine both reference trajactories to a dual reference trajectory
//
// const std::vector<TrajectoryStateOnSurface>& fwdTsosVec = fwdTraj->trajectoryStates();
// const std::vector<TrajectoryStateOnSurface>& bwdTsosVec = bwdTraj->trajectoryStates();
// theTsosVec.insert( theTsosVec.end(), fwdTsosVec.begin(), fwdTsosVec.end() );
// theTsosVec.insert( theTsosVec.end(), ++bwdTsosVec.begin(), bwdTsosVec.end() );
// Take hits as they come from the Kalman fit.
const ConstRecHitContainer &fwdRecHits = fwdTraj->recHits();
const ConstRecHitContainer &bwdRecHits = bwdTraj->recHits();
theRecHits.insert(theRecHits.end(), fwdRecHits.begin(), fwdRecHits.end());
theRecHits.insert(theRecHits.end(), ++bwdRecHits.begin(), bwdRecHits.end());
theParameters = this->extractParameters(refTsos);
unsigned int nParam = theParameters.num_row();
unsigned int nFwdMeas = nMeasPerHit * fwdTraj->numberOfHits();
unsigned int nBwdMeas = nMeasPerHit * bwdTraj->numberOfHits();
// theMeasurements.sub( 1, fwdTraj->measurements() );
// theMeasurements.sub( nFwdMeas+1, bwdTraj->measurements().sub( nMeasPerHit+1, nBwdMeas ) );
// theMeasurementsCov.sub( 1, fwdTraj->measurementErrors() );
// theMeasurementsCov.sub( nFwdMeas+1, bwdTraj->measurementErrors().sub( nMeasPerHit+1, nBwdMeas ) );
// theTrajectoryPositions.sub( 1, fwdTraj->trajectoryPositions() );
// theTrajectoryPositions.sub( nFwdMeas+1, bwdTraj->trajectoryPositions().sub( nMeasPerHit+1, nBwdMeas ) );
// theTrajectoryPositionCov.sub( 1, fwdTraj->trajectoryPositionErrors() );
// theTrajectoryPositionCov.sub( nFwdMeas+1, bwdTraj->trajectoryPositionErrors().sub( nMeasPerHit+1, nBwdMeas ) );
theDerivatives.sub(1, 1, fwdTraj->derivatives());
theDerivatives.sub(nFwdMeas+1, 1,
bwdTraj->derivatives().sub(nMeasPerHit+1, nBwdMeas, 1, nParam));
// FIXME: next lines stolen from ReferenceTrajectory, no idea whether OK or not...
if (refTsos.hasError()) {
AlgebraicSymMatrix parameterCov = asHepMatrix<5>(refTsos.localError().matrix());
theTrajectoryPositionCov = parameterCov.similarity(theDerivatives);
} else {
theTrajectoryPositionCov = AlgebraicSymMatrix(theDerivatives.num_row(), 1);
}
// Fill Kalman part, first for forward, then for backward part.
return this->fillKalmanPart(trajMeasurements, forwardRecHitNums, true, 0, residualMethod)
&& this->fillKalmanPart(trajMeasurements, backwardRecHitNums, false,
forwardRecHitNums.size(), residualMethod);
}
| ReferenceTrajectory * DualKalmanTrajectory::construct | ( | const Trajectory::DataContainer & | trajMeasurements, |
| const TrajectoryStateOnSurface & | referenceTsos, | ||
| const std::vector< unsigned int > & | recHits, | ||
| double | mass, | ||
| MaterialEffects | materialEffects, | ||
| const PropagationDirection | propDir, | ||
| const MagneticField * | magField, | ||
| bool | useBeamSpot, | ||
| const reco::BeamSpot & | beamSpot | ||
| ) | const [protected, virtual] |
Method to get a single ReferenceTrajectory for a half of the trajectory.
Definition at line 143 of file DualKalmanTrajectory.cc.
References i, and ReferenceTrajectoryBase::recHits().
{
ConstRecHitContainer recHits;
recHits.reserve(recHitNums.size());
for (unsigned int i = 0; i < recHitNums.size(); ++i) {
recHits.push_back(trajMeasurements[recHitNums[i]].recHit());
}
return new ReferenceTrajectory(referenceTsos, recHits, false, // hits are already ordered
magField, materialEffects, propDir, mass,
useBeamSpot, beamSpot);
}
| AlgebraicVector DualKalmanTrajectory::extractParameters | ( | const TrajectoryStateOnSurface & | referenceTsos | ) | const [protected, virtual] |
Definition at line 329 of file DualKalmanTrajectory.cc.
References TrajectoryStateOnSurface::localParameters(), and LocalTrajectoryParameters::mixedFormatVector().
Referenced by construct().
{
return asHepVector<5>( referenceTsos.localParameters().mixedFormatVector() );
}
| bool DualKalmanTrajectory::fillKalmanPart | ( | const Trajectory::DataContainer & | trajMeasurements, |
| const std::vector< unsigned int > & | recHitNums, | ||
| bool | startFirst, | ||
| unsigned int | iNextHit, | ||
| int | residualMethod | ||
| ) | [protected] |
Fill that part of data members that is different from DualReferenceTrajectory.
Definition at line 164 of file DualKalmanTrajectory.cc.
References Exception, fillMeasurementAndError1(), fillMeasurementAndError2(), fillTrajectoryPositions(), TrajectoryStateOnSurface::isValid(), TrajectoryMeasurement::recHit(), and ReferenceTrajectoryBase::theTsosVec.
Referenced by construct().
{
// startFirst==false: skip first hit of recHitNums
// iNextHit: first hit number to fill data members with
//
// Three approaches, choosen by 'residualMethod':
// 0: Just take hit error and the updated state as millepede needs it to redo the last fit step.
// 1: Use the unbiased residuals as for residual monitoring.
// 2: Use the _updated_ state and calculate the sigma that is part of
// the pull as sqrt(sigma_hit^2 - sigma_tsos^2).
// This should lead to the pull as defined on p. 236 of Blobel's book.
// Not sure whether this is 100% applicable/correct here, at least it might lead to
// numerical instabilities, cf.:
// https://hypernews.cern.ch/HyperNews/CMS/get/recoTracking/517/1/1/1/1/1.html
for (unsigned int iMeas = (startFirst ? 0 : 1); iMeas < recHitNums.size(); ++iMeas, ++iNextHit) {
const TrajectoryMeasurement &trajMeasurement = trajMeasurements[recHitNums[iMeas]];
TrajectoryStateOnSurface tsos; // depends on method
switch (residualMethod) {
case 0:
tsos = this->fillMeasurementAndError2(trajMeasurement.recHit(), iNextHit, trajMeasurement,
false); // plain hit error, not pull
break;
case 1:
tsos = this->fillMeasurementAndError1(trajMeasurement.recHit(), iNextHit, trajMeasurement);
break;
case 2:
tsos = this->fillMeasurementAndError2(trajMeasurement.recHit(), iNextHit, trajMeasurement,
true); // error of pull
break;
default:
throw cms::Exception("BadConfig")
<< "[DualKalmanTrajectory::fillKalmanPart] expect residualMethod == 0, 1 or 2, not "
<< residualMethod << ".";
}
if (!tsos.isValid()) return false;
theTsosVec.push_back(tsos);
this->fillTrajectoryPositions(trajMeasurement.recHit()->projectionMatrix(),
tsos, iNextHit);
}
return true;
}
| TrajectoryStateOnSurface DualKalmanTrajectory::fillMeasurementAndError1 | ( | const TransientTrackingRecHit::ConstRecHitPointer & | hitPtr, |
| unsigned int | iHit, | ||
| const TrajectoryMeasurement & | trajMeasurement | ||
| ) | [protected] |
helper for 'unbiased residual' method (i.e. returns merged fwd/bwd states)
Definition at line 216 of file DualKalmanTrajectory.cc.
References TrajectoryMeasurement::backwardPredictedState(), TrajectoryMeasurement::forwardPredictedState(), hitErrorWithAPE(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localError(), ReferenceTrajectoryBase::nMeasPerHit, LocalTrajectoryError::positionError(), ReferenceTrajectoryBase::theMeasurements, and ReferenceTrajectoryBase::theMeasurementsCov.
Referenced by fillKalmanPart().
{
// Get the measurements and their errors.
// We have to add error from hit and that tsos that is combination of fwd and bwd state.
TrajectoryStateCombiner tsosComb;
const TrajectoryStateOnSurface tsos = tsosComb(trajMeasurement.forwardPredictedState(),
trajMeasurement.backwardPredictedState());
if (tsos.isValid()) {
// No update of hit with tsos: it comes already from fwd+bwd tsos combination.
// See also https://hypernews.cern.ch/HyperNews/CMS/get/recoTracking/517/1.html .
// ConstRecHitPointer newHitPtr(hitPtr->canImproveWithTrack() ? hitPtr->clone(tsos) : hitPtr);
const LocalPoint localMeasurement(hitPtr->localPosition());
const LocalError hitErr(this->hitErrorWithAPE(hitPtr));
// tsos prediction includes APE of other hits:
const LocalError tsosErr(tsos.localError().positionError());
const LocalError localMeasurementCov(hitErr.xx() + tsosErr.xx(),
hitErr.xy() + tsosErr.xy(),
hitErr.yy() + tsosErr.yy());
theMeasurements[nMeasPerHit*iHit] = localMeasurement.x();
theMeasurements[nMeasPerHit*iHit+1] = localMeasurement.y();
theMeasurementsCov[nMeasPerHit*iHit] [nMeasPerHit*iHit] = localMeasurementCov.xx();
theMeasurementsCov[nMeasPerHit*iHit] [nMeasPerHit*iHit+1] = localMeasurementCov.xy();
theMeasurementsCov[nMeasPerHit*iHit+1][nMeasPerHit*iHit+1] = localMeasurementCov.yy();
}
return tsos;
}
| TrajectoryStateOnSurface DualKalmanTrajectory::fillMeasurementAndError2 | ( | const TransientTrackingRecHit::ConstRecHitPointer & | hitPtr, |
| unsigned int | iHit, | ||
| const TrajectoryMeasurement & | trajMeasurement, | ||
| bool | doPull | ||
| ) | [protected] |
helper for 'pull' (doPull=true) or 'hitError' method (doPull=false), returns updated tsos in both cases
Definition at line 251 of file DualKalmanTrajectory.cc.
References hitErrorWithAPE(), TrajectoryStateOnSurface::isValid(), TrajectoryStateOnSurface::localError(), ReferenceTrajectoryBase::nMeasPerHit, LocalTrajectoryError::positionError(), mathSSE::sqrt(), ReferenceTrajectoryBase::theMeasurements, ReferenceTrajectoryBase::theMeasurementsCov, and TrajectoryMeasurement::updatedState().
Referenced by fillKalmanPart().
{
// Get the measurements and their errors.
// We have to subtract error of updated tsos from hit.
const TrajectoryStateOnSurface tsos = trajMeasurement.updatedState();
if (tsos.isValid()) {
// No further update of hit:
// - The Kalman fit used hitPtr as it comes here (besides APE, see below).
// - If the hit errors improve, we might get (rare) problems of negative diagonal elements, see below.
// ConstRecHitPointer newHitPtr(hitPtr->canImproveWithTrack() ? hitPtr->clone(tsos) : hitPtr);
const LocalPoint localMeasurement(hitPtr->localPosition());
// Add APE manually to avoid that the hit error might be smaller than tsos error:
// - hit local errors are always without APE,
// - the tsos errors include APE since they come from track fit.
const LocalError hitErr(this->hitErrorWithAPE(hitPtr));
LocalError localMeasurementCov(hitErr.xx(), hitErr.xy(), hitErr.yy());
if (doPull) {
const LocalError tsosErr(tsos.localError().positionError());
// Should not be possible to become negative if all is correct - see above.
if (hitErr.xx() < tsosErr.xx() || hitErr.yy() < tsosErr.yy()) {
edm::LogError("Alignment") << "@SUB=DualKalmanTrajectory::fillMeasurementAndError2"
<< "not OK in subdet " << hitPtr->geographicalId().subdetId()
<< "\ns_x " << sqrt(hitErr.xx()) << " " << sqrt(tsosErr.xx())
<< "\ns_xy " << hitErr.xy() << " " << tsosErr.xy()
<< "\ns_y " << sqrt(hitErr.yy()) << " " << sqrt(tsosErr.yy());
return TrajectoryStateOnSurface(); // is invalid state
}
// cf. Blobel/Lohrmann, p. 236:
// But numerical stability? Should we return false if difference very small compared to values?
localMeasurementCov = LocalError(hitErr.xx() - tsosErr.xx(), // tsos puts correlation in,
hitErr.xy() - tsosErr.xy(), // even for 1D strip!
hitErr.yy() - tsosErr.yy());
}
theMeasurements[nMeasPerHit*iHit] = localMeasurement.x();
theMeasurements[nMeasPerHit*iHit+1] = localMeasurement.y();
theMeasurementsCov[nMeasPerHit*iHit] [nMeasPerHit*iHit] = localMeasurementCov.xx();
theMeasurementsCov[nMeasPerHit*iHit] [nMeasPerHit*iHit+1] = localMeasurementCov.xy();
theMeasurementsCov[nMeasPerHit*iHit+1][nMeasPerHit*iHit+1] = localMeasurementCov.yy();
}
return tsos;
}
| void DualKalmanTrajectory::fillTrajectoryPositions | ( | const AlgebraicMatrix & | projection, |
| const TrajectoryStateOnSurface & | tsos, | ||
| unsigned int | iHit | ||
| ) | [protected] |
fill trajectoryPositions
Definition at line 302 of file DualKalmanTrajectory.cc.
References TrajectoryStateOnSurface::localParameters(), LocalTrajectoryParameters::mixedFormatVector(), ReferenceTrajectoryBase::nMeasPerHit, and ReferenceTrajectoryBase::theTrajectoryPositions.
Referenced by fillKalmanPart().
{
// get the local coordinates of the reference trajectory
// (~copied from ReferenceTrajectory::fillTrajectoryPositions)
AlgebraicVector mixedLocalParams = asHepVector<5>(tsos.localParameters().mixedFormatVector());
const AlgebraicVector localPosition(projection * mixedLocalParams);
theTrajectoryPositions[nMeasPerHit*iHit] = localPosition[0];
theTrajectoryPositions[nMeasPerHit*iHit+1] = localPosition[1];
}
| LocalError DualKalmanTrajectory::hitErrorWithAPE | ( | const TransientTrackingRecHit::ConstRecHitPointer & | hitPtr | ) | const [protected] |
local error including APE if APE is on
Definition at line 316 of file DualKalmanTrajectory.cc.
References alignCSCRings::help, and HelpertRecHit2DLocalPos::parError().
Referenced by fillMeasurementAndError1(), and fillMeasurementAndError2().
{
if (hitPtr->det() && hitPtr->det()->alignmentPositionError()) {
HelpertRecHit2DLocalPos help;
const AlgebraicSymMatrix errMat(help.parError(hitPtr->localPositionError(), *(hitPtr->det())));
return LocalError(errMat[0][0], errMat[0][1], errMat[1][1]);
} else {
return hitPtr->localPositionError();
}
}
| const PropagationDirection DualKalmanTrajectory::oppositeDirection | ( | const PropagationDirection | propDir | ) | const [inline, protected] |
Definition at line 102 of file DualKalmanTrajectory.h.
References alongMomentum, anyDirection, and oppositeToMomentum.
Referenced by construct().
{ return ((propDir == anyDirection) ? anyDirection
: ((propDir == alongMomentum) ? oppositeToMomentum : alongMomentum));
}
1.7.3