TwoBodyDecayEstimator Class Reference

#include <TwoBodyDecayEstimator.h>

Public Types
typedef PerigeeLinearizedTrackState::RefCountedLinearizedTrackState	RefCountedLinearizedTrackState

Public Member Functions
virtual TwoBodyDecayEstimator *	clone (void) const
virtual TwoBodyDecay	estimate (const std::vector< RefCountedLinearizedTrackState > &linTracks, const TwoBodyDecayParameters &linearizationPoint, const TwoBodyDecayVirtualMeasurement &vm) const
int	ndf (void) const
const AlgebraicVector &	pulls (void) const
	TwoBodyDecayEstimator (const edm::ParameterSet &config)
virtual	~TwoBodyDecayEstimator (void)

Protected Member Functions
virtual bool	constructMatrices (const std::vector< RefCountedLinearizedTrackState > &linTracks, const TwoBodyDecayParameters &linearizationPoint, const TwoBodyDecayVirtualMeasurement &vm, AlgebraicVector &vecM, AlgebraicSymMatrix &matG, AlgebraicMatrix &matA) const

Private Member Functions
bool	checkValues (const AlgebraicVector &vec) const

Private Attributes
int	theMaxIterations
double	theMaxIterDiff
int	theNdf
AlgebraicVector	thePulls
double	theRobustificationConstant
bool	theUseInvariantMass

Detailed Description

Class TwoBodyDecayEstimator estimates the decay parameters and the corresponding error matrix (see TwoBodyDecayParameter) from two linearized tracks, and an hypothesis for the primary particle's mass and width as well as the mass of the secondary particles. It utilizes a robust M-estimator.

/author Edmund Widl

Definition at line 20 of file TwoBodyDecayEstimator.h.

Member Typedef Documentation

RefCountedLinearizedTrackState

typedef PerigeeLinearizedTrackState::RefCountedLinearizedTrackState TwoBodyDecayEstimator::RefCountedLinearizedTrackState

Definition at line 22 of file TwoBodyDecayEstimator.h.

Constructor & Destructor Documentation

TwoBodyDecayEstimator()

TwoBodyDecayEstimator::TwoBodyDecayEstimator

(

const edm::ParameterSet &

config

)

Definition at line 11 of file TwoBodyDecayEstimator.cc.

                                                                          : theNdf(0) {
  const edm::ParameterSet &estimatorConfig = config.getParameter<edm::ParameterSet>("EstimatorParameters");
 
  theRobustificationConstant = estimatorConfig.getUntrackedParameter<double>("RobustificationConstant", 1.0);
  theMaxIterDiff = estimatorConfig.getUntrackedParameter<double>("MaxIterationDifference", 1e-2);
  theMaxIterations = estimatorConfig.getUntrackedParameter<int>("MaxIterations", 100);
  theUseInvariantMass = estimatorConfig.getUntrackedParameter<bool>("UseInvariantMass", true);
}

References MillePedeFileConverter_cfg::e, edm::ParameterSet::getUntrackedParameter(), theMaxIterations, theMaxIterDiff, theRobustificationConstant, and theUseInvariantMass.

Referenced by clone().

~TwoBodyDecayEstimator()

virtual TwoBodyDecayEstimator::~TwoBodyDecayEstimator ( void  )

inlinevirtual

Definition at line 25 of file TwoBodyDecayEstimator.h.

{}

Member Function Documentation

checkValues()

bool TwoBodyDecayEstimator::checkValues ( const AlgebraicVector &  vec ) const

private

Definition at line 215 of file TwoBodyDecayEstimator.cc.

                                                                        {
  bool isNotFinite = false;
 
  for (int i = 0; i < vec.num_col(); ++i)
    isNotFinite |= edm::isNotFinite(vec[i]);
 
  return isNotFinite;
}

References mps_fire::i, and edm::isNotFinite().

Referenced by constructMatrices().

clone()

virtual TwoBodyDecayEstimator* TwoBodyDecayEstimator::clone ( void  ) const

inlinevirtual

Definition at line 34 of file TwoBodyDecayEstimator.h.

{ return new TwoBodyDecayEstimator(*this); }

References TwoBodyDecayEstimator().

constructMatrices()

bool TwoBodyDecayEstimator::constructMatrices ( const std::vector< RefCountedLinearizedTrackState > &  linTracks, const TwoBodyDecayParameters &  linearizationPoint, const TwoBodyDecayVirtualMeasurement &  vm, AlgebraicVector &  vecM, AlgebraicSymMatrix &  matG, AlgebraicMatrix &  matA  ) const

protectedvirtual

Definition at line 99 of file TwoBodyDecayEstimator.cc.

                                                                           {
  PerigeeLinearizedTrackState *linTrack1 = dynamic_cast<PerigeeLinearizedTrackState *>(linTracks[0].get());
  PerigeeLinearizedTrackState *linTrack2 = dynamic_cast<PerigeeLinearizedTrackState *>(linTracks[1].get());
 
  if (!linTrack1 || !linTrack2)
    return false;
 
  AlgebraicVector trackParam1 = asHepVector(linTrack1->predictedStateParameters());
  AlgebraicVector trackParam2 = asHepVector(linTrack2->predictedStateParameters());
 
  if (checkValues(trackParam1) || checkValues(trackParam2) || checkValues(linearizationPoint.parameters()))
    return false;
 
  AlgebraicVector vecLinParam = linearizationPoint.sub(TwoBodyDecayParameters::px, TwoBodyDecayParameters::mass);
 
  double zMagField = linTrack1->track().field()->inInverseGeV(linTrack1->linearizationPoint()).z();
 
  int checkInversion = 0;
 
  TwoBodyDecayDerivatives tpeDerivatives(linearizationPoint[TwoBodyDecayParameters::mass], vm.secondaryMass());
  std::pair<AlgebraicMatrix, AlgebraicMatrix> derivatives = tpeDerivatives.derivatives(linearizationPoint);
 
  TwoBodyDecayModel decayModel(linearizationPoint[TwoBodyDecayParameters::mass], vm.secondaryMass());
  std::pair<AlgebraicVector, AlgebraicVector> linCartMomenta = decayModel.cartesianSecondaryMomenta(linearizationPoint);
 
  // first track
  AlgebraicMatrix matA1 = asHepMatrix(linTrack1->positionJacobian());
  AlgebraicMatrix matB1 = asHepMatrix(linTrack1->momentumJacobian());
  AlgebraicVector vecC1 = asHepVector(linTrack1->constantTerm());
 
  AlgebraicVector curvMomentum1 = asHepVector(linTrack1->predictedStateMomentumParameters());
  AlgebraicMatrix curv2cart1 = decayModel.curvilinearToCartesianJacobian(curvMomentum1, zMagField);
 
  AlgebraicVector cartMomentum1 = decayModel.convertCurvilinearToCartesian(curvMomentum1, zMagField);
  vecC1 += matB1 * (curvMomentum1 - curv2cart1 * cartMomentum1);
  matB1 = matB1 * curv2cart1;
 
  AlgebraicMatrix matF1 = derivatives.first;
  AlgebraicVector vecD1 = linCartMomenta.first - matF1 * vecLinParam;
  AlgebraicVector vecM1 = trackParam1 - vecC1 - matB1 * vecD1;
  AlgebraicSymMatrix covM1 = asHepMatrix(linTrack1->predictedStateError());
 
  AlgebraicSymMatrix matG1 = covM1.inverse(checkInversion);
  if (checkInversion != 0) {
    LogDebug("Alignment") << "@SUB=TwoBodyDecayEstimator::constructMatrices"
                          << "Matrix covM1 not invertible.";
    return false;
  }
 
  // diagonalize for robustification
  AlgebraicMatrix matU1 = diagonalize(&matG1).T();
 
  // second track
  AlgebraicMatrix matA2 = asHepMatrix(linTrack2->positionJacobian());
  AlgebraicMatrix matB2 = asHepMatrix(linTrack2->momentumJacobian());
  AlgebraicVector vecC2 = asHepVector(linTrack2->constantTerm());
 
  AlgebraicVector curvMomentum2 = asHepVector(linTrack2->predictedStateMomentumParameters());
  AlgebraicMatrix curv2cart2 = decayModel.curvilinearToCartesianJacobian(curvMomentum2, zMagField);
 
  AlgebraicVector cartMomentum2 = decayModel.convertCurvilinearToCartesian(curvMomentum2, zMagField);
  vecC2 += matB2 * (curvMomentum2 - curv2cart2 * cartMomentum2);
  matB2 = matB2 * curv2cart2;
 
  AlgebraicMatrix matF2 = derivatives.second;
  AlgebraicVector vecD2 = linCartMomenta.second - matF2 * vecLinParam;
  AlgebraicVector vecM2 = trackParam2 - vecC2 - matB2 * vecD2;
  AlgebraicSymMatrix covM2 = asHepMatrix(linTrack2->predictedStateError());
 
  AlgebraicSymMatrix matG2 = covM2.inverse(checkInversion);
  if (checkInversion != 0) {
    LogDebug("Alignment") << "@SUB=TwoBodyDecayEstimator::constructMatrices"
                          << "Matrix covM2 not invertible.";
    return false;
  }
 
  // diagonalize for robustification
  AlgebraicMatrix matU2 = diagonalize(&matG2).T();
 
  // combine first and second track
  vecM = AlgebraicVector(14, 0);
  vecM.sub(1, matU1 * vecM1);
  vecM.sub(6, matU2 * vecM2);
  // virtual measurement of the primary mass
  vecM(11) = vm.primaryMass();
  // virtual measurement of the beam spot
  vecM.sub(12, vm.beamSpotPosition());
 
  // full weight matrix
  matG = AlgebraicSymMatrix(14, 0);
  matG.sub(1, matG1);
  matG.sub(6, matG2);
  // virtual measurement error of the primary mass
  matG[10][10] = 1. / (vm.primaryWidth() * vm.primaryWidth());
  // virtual measurement error of the beam spot
  matG.sub(12, vm.beamSpotError().inverse(checkInversion));
 
  // full derivative matrix
  matA = AlgebraicMatrix(14, 9, 0);
  matA.sub(1, 1, matU1 * matA1);
  matA.sub(6, 1, matU2 * matA2);
  matA.sub(1, 4, matU1 * matB1 * matF1);
  matA.sub(6, 4, matU2 * matB2 * matF2);
  matA(11, 9) = 1.;  // mass
  matA(12, 1) = 1.;  // vx
  matA(13, 2) = 1.;  // vy
  matA(14, 3) = 1.;  // vz
 
  return true;
}

References asHepMatrix(), asHepVector(), TwoBodyDecayVirtualMeasurement::beamSpotError(), TwoBodyDecayVirtualMeasurement::beamSpotPosition(), TwoBodyDecayModel::cartesianSecondaryMomenta(), checkValues(), PerigeeLinearizedTrackState::constantTerm(), TwoBodyDecayDerivatives::derivatives(), reco::TransientTrack::field(), get, MagneticField::inInverseGeV(), PerigeeLinearizedTrackState::linearizationPoint(), LogDebug, TwoBodyDecayParameters::mass, PerigeeLinearizedTrackState::momentumJacobian(), TwoBodyDecayParameters::parameters(), PerigeeLinearizedTrackState::positionJacobian(), PerigeeLinearizedTrackState::predictedStateError(), PerigeeLinearizedTrackState::predictedStateMomentumParameters(), PerigeeLinearizedTrackState::predictedStateParameters(), TwoBodyDecayVirtualMeasurement::primaryMass(), TwoBodyDecayVirtualMeasurement::primaryWidth(), TwoBodyDecayParameters::px, TwoBodyDecayVirtualMeasurement::secondaryMass(), TwoBodyDecayParameters::sub(), PerigeeLinearizedTrackState::track(), and z.

Referenced by estimate().

estimate()

TwoBodyDecay TwoBodyDecayEstimator::estimate ( const std::vector< RefCountedLinearizedTrackState > &  linTracks, const TwoBodyDecayParameters &  linearizationPoint, const TwoBodyDecayVirtualMeasurement &  vm  ) const

virtual

Definition at line 20 of file TwoBodyDecayEstimator.cc.

                                                                                             {
  if (linTracks.size() != 2) {
    edm::LogInfo("Alignment") << "@SUB=TwoBodyDecayEstimator::estimate"
                              << "Need 2 linearized tracks, got " << linTracks.size() << ".\n";
    return TwoBodyDecay();
  }
 
  AlgebraicVector vecM;
  AlgebraicSymMatrix matG;
  AlgebraicMatrix matA;
 
  bool check = constructMatrices(linTracks, linearizationPoint, vm, vecM, matG, matA);
  if (!check)
    return TwoBodyDecay();
 
  AlgebraicSymMatrix matGPrime;
  AlgebraicSymMatrix invAtGPrimeA;
  AlgebraicVector vecEstimate;
  AlgebraicVector res;
 
  int nIterations = 0;
  bool stopIteration = false;
 
  // initialization - values are never used
  int checkInversion = 0;
  double chi2 = 0.;
  double oldChi2 = 0.;
  bool isValid = true;
 
  while (!stopIteration) {
    matGPrime = matG;
 
    // compute weights
    if (nIterations > 0) {
      for (int i = 0; i < 10; i++) {
        double sigma = 1. / sqrt(matG[i][i]);
        double sigmaTimesR = sigma * theRobustificationConstant;
        double absRes = fabs(res[i]);
        if (absRes > sigmaTimesR)
          matGPrime[i][i] *= sigmaTimesR / absRes;
      }
    }
 
    // make LS-fit
    invAtGPrimeA = (matGPrime.similarityT(matA)).inverse(checkInversion);
    if (checkInversion != 0) {
      LogDebug("Alignment") << "@SUB=TwoBodyDecayEstimator::estimate"
                            << "Matrix At*G'*A not invertible (in iteration " << nIterations
                            << ", ifail = " << checkInversion << ").\n";
      isValid = false;
      break;
    }
    vecEstimate = invAtGPrimeA * matA.T() * matGPrime * vecM;
    res = matA * vecEstimate - vecM;
    chi2 = dot(res, matGPrime * res);
 
    if ((nIterations > 0) && (fabs(chi2 - oldChi2) < theMaxIterDiff))
      stopIteration = true;
    if (nIterations == theMaxIterations)
      stopIteration = true;
 
    oldChi2 = chi2;
    nIterations++;
  }
 
  if (isValid) {
    AlgebraicSymMatrix pullsCov = matGPrime.inverse(checkInversion) - invAtGPrimeA.similarity(matA);
    thePulls = AlgebraicVector(matG.num_col(), 0);
    for (int i = 0; i < pullsCov.num_col(); i++)
      thePulls[i] = res[i] / sqrt(pullsCov[i][i]);
  }
 
  theNdf = matA.num_row() - matA.num_col();
 
  return TwoBodyDecay(TwoBodyDecayParameters(vecEstimate, invAtGPrimeA), chi2, isValid, vm);
}

References RPCNoise_example::check, hltPixelTracks_cff::chi2, constructMatrices(), dot(), mps_fire::i, LogDebug, mathSSE::sqrt(), theMaxIterations, theMaxIterDiff, theNdf, thePulls, and theRobustificationConstant.

Referenced by TwoBodyDecayFitter::estimate().

ndf()

int TwoBodyDecayEstimator::ndf ( void  ) const

inline

Definition at line 31 of file TwoBodyDecayEstimator.h.

{ return theNdf; }

References theNdf.

pulls()

const AlgebraicVector& TwoBodyDecayEstimator::pulls ( void  ) const

inline

Definition at line 32 of file TwoBodyDecayEstimator.h.

{ return thePulls; }

References thePulls.

Member Data Documentation

theMaxIterations

int TwoBodyDecayEstimator::theMaxIterations

private

Definition at line 49 of file TwoBodyDecayEstimator.h.

Referenced by estimate(), and TwoBodyDecayEstimator().

theMaxIterDiff

double TwoBodyDecayEstimator::theMaxIterDiff

private

Definition at line 48 of file TwoBodyDecayEstimator.h.

Referenced by estimate(), and TwoBodyDecayEstimator().

theNdf

int TwoBodyDecayEstimator::theNdf

mutableprivate

Definition at line 52 of file TwoBodyDecayEstimator.h.

Referenced by estimate(), and ndf().

thePulls

AlgebraicVector TwoBodyDecayEstimator::thePulls

mutableprivate

Definition at line 53 of file TwoBodyDecayEstimator.h.

Referenced by estimate(), and pulls().

theRobustificationConstant

double TwoBodyDecayEstimator::theRobustificationConstant

private

Definition at line 47 of file TwoBodyDecayEstimator.h.

Referenced by estimate(), and TwoBodyDecayEstimator().

theUseInvariantMass

bool TwoBodyDecayEstimator::theUseInvariantMass

private

Definition at line 50 of file TwoBodyDecayEstimator.h.

Referenced by TwoBodyDecayEstimator().