CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes
TwoBodyDecayEstimator Class Reference

#include <TwoBodyDecayEstimator.h>

Public Types

typedef PerigeeLinearizedTrackState::RefCountedLinearizedTrackState RefCountedLinearizedTrackState
 

Public Member Functions

virtual TwoBodyDecayEstimatorclone (void) const
 
virtual TwoBodyDecay estimate (const std::vector< RefCountedLinearizedTrackState > &linTracks, const TwoBodyDecayParameters &linearizationPoint, const TwoBodyDecayVirtualMeasurement &vm) const
 
int ndf (void) const
 
const AlgebraicVectorpulls (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

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.

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

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.

25 {}

Member Function Documentation

◆ checkValues()

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

Definition at line 215 of file TwoBodyDecayEstimator.cc.

215  {
216  bool isNotFinite = false;
217 
218  for (int i = 0; i < vec.num_col(); ++i)
219  isNotFinite |= edm::isNotFinite(vec[i]);
220 
221  return isNotFinite;
222 }

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.

34 { 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.

104  {
105  PerigeeLinearizedTrackState *linTrack1 = dynamic_cast<PerigeeLinearizedTrackState *>(linTracks[0].get());
106  PerigeeLinearizedTrackState *linTrack2 = dynamic_cast<PerigeeLinearizedTrackState *>(linTracks[1].get());
107 
108  if (!linTrack1 || !linTrack2)
109  return false;
110 
111  AlgebraicVector trackParam1 = asHepVector(linTrack1->predictedStateParameters());
112  AlgebraicVector trackParam2 = asHepVector(linTrack2->predictedStateParameters());
113 
114  if (checkValues(trackParam1) || checkValues(trackParam2) || checkValues(linearizationPoint.parameters()))
115  return false;
116 
118 
119  double zMagField = linTrack1->track().field()->inInverseGeV(linTrack1->linearizationPoint()).z();
120 
121  int checkInversion = 0;
122 
123  TwoBodyDecayDerivatives tpeDerivatives(linearizationPoint[TwoBodyDecayParameters::mass], vm.secondaryMass());
124  std::pair<AlgebraicMatrix, AlgebraicMatrix> derivatives = tpeDerivatives.derivatives(linearizationPoint);
125 
126  TwoBodyDecayModel decayModel(linearizationPoint[TwoBodyDecayParameters::mass], vm.secondaryMass());
127  std::pair<AlgebraicVector, AlgebraicVector> linCartMomenta = decayModel.cartesianSecondaryMomenta(linearizationPoint);
128 
129  // first track
130  AlgebraicMatrix matA1 = asHepMatrix(linTrack1->positionJacobian());
131  AlgebraicMatrix matB1 = asHepMatrix(linTrack1->momentumJacobian());
132  AlgebraicVector vecC1 = asHepVector(linTrack1->constantTerm());
133 
134  AlgebraicVector curvMomentum1 = asHepVector(linTrack1->predictedStateMomentumParameters());
135  AlgebraicMatrix curv2cart1 = decayModel.curvilinearToCartesianJacobian(curvMomentum1, zMagField);
136 
137  AlgebraicVector cartMomentum1 = decayModel.convertCurvilinearToCartesian(curvMomentum1, zMagField);
138  vecC1 += matB1 * (curvMomentum1 - curv2cart1 * cartMomentum1);
139  matB1 = matB1 * curv2cart1;
140 
141  AlgebraicMatrix matF1 = derivatives.first;
142  AlgebraicVector vecD1 = linCartMomenta.first - matF1 * vecLinParam;
143  AlgebraicVector vecM1 = trackParam1 - vecC1 - matB1 * vecD1;
144  AlgebraicSymMatrix covM1 = asHepMatrix(linTrack1->predictedStateError());
145 
146  AlgebraicSymMatrix matG1 = covM1.inverse(checkInversion);
147  if (checkInversion != 0) {
148  LogDebug("Alignment") << "@SUB=TwoBodyDecayEstimator::constructMatrices"
149  << "Matrix covM1 not invertible.";
150  return false;
151  }
152 
153  // diagonalize for robustification
154  AlgebraicMatrix matU1 = diagonalize(&matG1).T();
155 
156  // second track
157  AlgebraicMatrix matA2 = asHepMatrix(linTrack2->positionJacobian());
158  AlgebraicMatrix matB2 = asHepMatrix(linTrack2->momentumJacobian());
159  AlgebraicVector vecC2 = asHepVector(linTrack2->constantTerm());
160 
161  AlgebraicVector curvMomentum2 = asHepVector(linTrack2->predictedStateMomentumParameters());
162  AlgebraicMatrix curv2cart2 = decayModel.curvilinearToCartesianJacobian(curvMomentum2, zMagField);
163 
164  AlgebraicVector cartMomentum2 = decayModel.convertCurvilinearToCartesian(curvMomentum2, zMagField);
165  vecC2 += matB2 * (curvMomentum2 - curv2cart2 * cartMomentum2);
166  matB2 = matB2 * curv2cart2;
167 
168  AlgebraicMatrix matF2 = derivatives.second;
169  AlgebraicVector vecD2 = linCartMomenta.second - matF2 * vecLinParam;
170  AlgebraicVector vecM2 = trackParam2 - vecC2 - matB2 * vecD2;
171  AlgebraicSymMatrix covM2 = asHepMatrix(linTrack2->predictedStateError());
172 
173  AlgebraicSymMatrix matG2 = covM2.inverse(checkInversion);
174  if (checkInversion != 0) {
175  LogDebug("Alignment") << "@SUB=TwoBodyDecayEstimator::constructMatrices"
176  << "Matrix covM2 not invertible.";
177  return false;
178  }
179 
180  // diagonalize for robustification
181  AlgebraicMatrix matU2 = diagonalize(&matG2).T();
182 
183  // combine first and second track
184  vecM = AlgebraicVector(14, 0);
185  vecM.sub(1, matU1 * vecM1);
186  vecM.sub(6, matU2 * vecM2);
187  // virtual measurement of the primary mass
188  vecM(11) = vm.primaryMass();
189  // virtual measurement of the beam spot
190  vecM.sub(12, vm.beamSpotPosition());
191 
192  // full weight matrix
193  matG = AlgebraicSymMatrix(14, 0);
194  matG.sub(1, matG1);
195  matG.sub(6, matG2);
196  // virtual measurement error of the primary mass
197  matG[10][10] = 1. / (vm.primaryWidth() * vm.primaryWidth());
198  // virtual measurement error of the beam spot
199  matG.sub(12, vm.beamSpotError().inverse(checkInversion));
200 
201  // full derivative matrix
202  matA = AlgebraicMatrix(14, 9, 0);
203  matA.sub(1, 1, matU1 * matA1);
204  matA.sub(6, 1, matU2 * matA2);
205  matA.sub(1, 4, matU1 * matB1 * matF1);
206  matA.sub(6, 4, matU2 * matB2 * matF2);
207  matA(11, 9) = 1.; // mass
208  matA(12, 1) = 1.; // vx
209  matA(13, 2) = 1.; // vy
210  matA(14, 3) = 1.; // vz
211 
212  return true;
213 }

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.

22  {
23  if (linTracks.size() != 2) {
24  edm::LogInfo("Alignment") << "@SUB=TwoBodyDecayEstimator::estimate"
25  << "Need 2 linearized tracks, got " << linTracks.size() << ".\n";
26  return TwoBodyDecay();
27  }
28 
29  AlgebraicVector vecM;
30  AlgebraicSymMatrix matG;
31  AlgebraicMatrix matA;
32 
33  bool check = constructMatrices(linTracks, linearizationPoint, vm, vecM, matG, matA);
34  if (!check)
35  return TwoBodyDecay();
36 
37  AlgebraicSymMatrix matGPrime;
38  AlgebraicSymMatrix invAtGPrimeA;
39  AlgebraicVector vecEstimate;
41 
42  int nIterations = 0;
43  bool stopIteration = false;
44 
45  // initialization - values are never used
46  int checkInversion = 0;
47  double chi2 = 0.;
48  double oldChi2 = 0.;
49  bool isValid = true;
50 
51  while (!stopIteration) {
52  matGPrime = matG;
53 
54  // compute weights
55  if (nIterations > 0) {
56  for (int i = 0; i < 10; i++) {
57  double sigma = 1. / sqrt(matG[i][i]);
58  double sigmaTimesR = sigma * theRobustificationConstant;
59  double absRes = fabs(res[i]);
60  if (absRes > sigmaTimesR)
61  matGPrime[i][i] *= sigmaTimesR / absRes;
62  }
63  }
64 
65  // make LS-fit
66  invAtGPrimeA = (matGPrime.similarityT(matA)).inverse(checkInversion);
67  if (checkInversion != 0) {
68  LogDebug("Alignment") << "@SUB=TwoBodyDecayEstimator::estimate"
69  << "Matrix At*G'*A not invertible (in iteration " << nIterations
70  << ", ifail = " << checkInversion << ").\n";
71  isValid = false;
72  break;
73  }
74  vecEstimate = invAtGPrimeA * matA.T() * matGPrime * vecM;
75  res = matA * vecEstimate - vecM;
76  chi2 = dot(res, matGPrime * res);
77 
78  if ((nIterations > 0) && (fabs(chi2 - oldChi2) < theMaxIterDiff))
79  stopIteration = true;
80  if (nIterations == theMaxIterations)
81  stopIteration = true;
82 
83  oldChi2 = chi2;
84  nIterations++;
85  }
86 
87  if (isValid) {
88  AlgebraicSymMatrix pullsCov = matGPrime.inverse(checkInversion) - invAtGPrimeA.similarity(matA);
89  thePulls = AlgebraicVector(matG.num_col(), 0);
90  for (int i = 0; i < pullsCov.num_col(); i++)
91  thePulls[i] = res[i] / sqrt(pullsCov[i][i]);
92  }
93 
94  theNdf = matA.num_row() - matA.num_col();
95 
96  return TwoBodyDecay(TwoBodyDecayParameters(vecEstimate, invAtGPrimeA), chi2, isValid, vm);
97 }

References RPCNoise_example::check, hltPixelTracks_cff::chi2, constructMatrices(), dot(), mps_fire::i, sistrip::SpyUtilities::isValid(), 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.

31 { return theNdf; }

References theNdf.

◆ pulls()

const AlgebraicVector& TwoBodyDecayEstimator::pulls ( void  ) const
inline

Definition at line 32 of file TwoBodyDecayEstimator.h.

32 { 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().

dot
T dot(const Basic3DVector &v) const
Scalar product, or "dot" product, with a vector of same type.
Definition: Basic3DVectorLD.h:212
PerigeeLinearizedTrackState::linearizationPoint
const GlobalPoint & linearizationPoint() const override
Definition: PerigeeLinearizedTrackState.h:50
TwoBodyDecayModel::cartesianSecondaryMomenta
const std::pair< AlgebraicVector, AlgebraicVector > cartesianSecondaryMomenta(const AlgebraicVector &param)
Definition: TwoBodyDecayModel.cc:80
mps_fire.i
i
Definition: mps_fire.py:428
PerigeeLinearizedTrackState::track
reco::TransientTrack track() const override
Definition: PerigeeLinearizedTrackState.h:52
TwoBodyDecayEstimator::checkValues
bool checkValues(const AlgebraicVector &vec) const
Definition: TwoBodyDecayEstimator.cc:215
edm::isNotFinite
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
TwoBodyDecayParameters::sub
const AlgebraicVector sub(ParameterName first, ParameterName last) const
Get specified range of decay parameters.
Definition: TwoBodyDecayParameters.h:46
TwoBodyDecayParameters::mass
Definition: TwoBodyDecayParameters.h:17
TwoBodyDecayDerivatives::derivatives
const std::pair< AlgebraicMatrix, AlgebraicMatrix > derivatives(const TwoBodyDecay &tbd) const
Definition: TwoBodyDecayDerivatives.cc:14
TwoBodyDecayEstimator::theMaxIterDiff
double theMaxIterDiff
Definition: TwoBodyDecayEstimator.h:48
TwoBodyDecayParameters::parameters
const AlgebraicVector & parameters(void) const
Get decay parameters.
Definition: TwoBodyDecayParameters.h:34
PerigeeLinearizedTrackState::constantTerm
const AlgebraicVector5 & constantTerm() const override
Definition: PerigeeLinearizedTrackState.h:176
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
hltPixelTracks_cff.chi2
chi2
Definition: hltPixelTracks_cff.py:25
PerigeeLinearizedTrackState::momentumJacobian
const AlgebraicMatrix53 & momentumJacobian() const override
Definition: PerigeeLinearizedTrackState.h:194
TwoBodyDecayVirtualMeasurement::primaryWidth
const double & primaryWidth(void) const
Definition: TwoBodyDecayVirtualMeasurement.h:37
PerigeeLinearizedTrackState::positionJacobian
const AlgebraicMatrix53 & positionJacobian() const override
Definition: PerigeeLinearizedTrackState.h:185
config
Definition: config.py:1
RPCNoise_example.check
check
Definition: RPCNoise_example.py:71
AlgebraicVector
CLHEP::HepVector AlgebraicVector
Definition: AlgebraicObjects.h:13
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
TwoBodyDecayEstimator::TwoBodyDecayEstimator
TwoBodyDecayEstimator(const edm::ParameterSet &config)
Definition: TwoBodyDecayEstimator.cc:11
DDAxes::z
TwoBodyDecayEstimator::theMaxIterations
int theMaxIterations
Definition: TwoBodyDecayEstimator.h:49
TwoBodyDecayEstimator::theRobustificationConstant
double theRobustificationConstant
Definition: TwoBodyDecayEstimator.h:47
sistrip::SpyUtilities::isValid
const bool isValid(const Frame &aFrame, const FrameQuality &aQuality, const uint16_t aExpectedPos)
Definition: SiStripSpyUtilities.cc:124
MagneticField::inInverseGeV
GlobalVector inInverseGeV(const GlobalPoint &gp) const
Field value ad specified global point, in 1/Gev.
Definition: MagneticField.h:36
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
edm::ParameterSet
Definition: ParameterSet.h:47
TwoBodyDecayDerivatives
Definition: TwoBodyDecayDerivatives.h:14
TwoBodyDecayModel
Definition: TwoBodyDecayModel.h:15
TwoBodyDecayVirtualMeasurement::beamSpotPosition
const AlgebraicVector beamSpotPosition(void) const
Definition: TwoBodyDecayVirtualMeasurement.h:41
PerigeeLinearizedTrackState::predictedStateMomentumParameters
AlgebraicVector3 predictedStateMomentumParameters() const override
Definition: PerigeeLinearizedTrackState.h:230
TwoBodyDecayParameters
Definition: TwoBodyDecayParameters.h:14
AlgebraicSymMatrix
CLHEP::HepSymMatrix AlgebraicSymMatrix
Definition: AlgebraicObjects.h:15
get
#define get
res
Definition: Electron.h:6
reco::TransientTrack::field
const MagneticField * field() const
Definition: TransientTrack.h:108
TwoBodyDecayVirtualMeasurement::beamSpotError
const AlgebraicSymMatrix beamSpotError(void) const
Definition: TwoBodyDecayVirtualMeasurement.h:42
TwoBodyDecayVirtualMeasurement::primaryMass
const double & primaryMass(void) const
Definition: TwoBodyDecayVirtualMeasurement.h:36
TwoBodyDecayVirtualMeasurement::secondaryMass
const double & secondaryMass(void) const
Definition: TwoBodyDecayVirtualMeasurement.h:38
AlgebraicMatrix
CLHEP::HepMatrix AlgebraicMatrix
Definition: AlgebraicObjects.h:14
TwoBodyDecayParameters::px
Definition: TwoBodyDecayParameters.h:17
TwoBodyDecayEstimator::theUseInvariantMass
bool theUseInvariantMass
Definition: TwoBodyDecayEstimator.h:50
TwoBodyDecay
Definition: TwoBodyDecay.h:15
asHepVector
CLHEP::HepVector asHepVector(const ROOT::Math::SVector< double, N > &v)
Definition: Migration.h:53
TwoBodyDecayEstimator::theNdf
int theNdf
Definition: TwoBodyDecayEstimator.h:52
TwoBodyDecayEstimator::thePulls
AlgebraicVector thePulls
Definition: TwoBodyDecayEstimator.h:53
TwoBodyDecayEstimator::constructMatrices
virtual bool constructMatrices(const std::vector< RefCountedLinearizedTrackState > &linTracks, const TwoBodyDecayParameters &linearizationPoint, const TwoBodyDecayVirtualMeasurement &vm, AlgebraicVector &vecM, AlgebraicSymMatrix &matG, AlgebraicMatrix &matA) const
Definition: TwoBodyDecayEstimator.cc:99
PerigeeLinearizedTrackState
Definition: PerigeeLinearizedTrackState.h:33
PerigeeLinearizedTrackState::predictedStateError
AlgebraicSymMatrix55 predictedStateError() const override
Definition: PerigeeLinearizedTrackState.h:247
asHepMatrix
CLHEP::HepMatrix asHepMatrix(const ROOT::Math::SMatrix< double, N1, N2, typename ROOT::Math::MatRepStd< double, N1, N2 > > &rm)
Definition: Migration.h:60
PerigeeLinearizedTrackState::predictedStateParameters
AlgebraicVector5 predictedStateParameters() const override
Definition: PerigeeLinearizedTrackState.h:224
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37