CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
List of all members | Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes
ConstrainedTreeBuilderT Class Reference

#include <ConstrainedTreeBuilderT.h>

Public Member Functions

template<int nTrk>
RefCountedKinematicTree buildTree (const std::vector< RefCountedKinematicParticle > &initialParticles, const std::vector< KinematicState > &finalStates, const RefCountedKinematicVertex vtx, const ROOT::Math::SMatrix< double, 3+7 *nTrk, 3+7 *nTrk, ROOT::Math::MatRepSym< double, 3+7 *nTrk > > &fCov) const
 
 ConstrainedTreeBuilderT ()
 
 ~ConstrainedTreeBuilderT ()
 

Private Member Functions

RefCountedKinematicTree buildRealTree (const RefCountedKinematicParticle virtualParticle, const RefCountedKinematicVertex vtx, const std::vector< RefCountedKinematicParticle > &particles) const
 

Static Private Member Functions

template<int nTrk>
static AlgebraicSymMatrix77 covarianceMatrix (const std::vector< RefCountedKinematicParticle > &rPart, const AlgebraicVector7 &newPar, const ROOT::Math::SMatrix< double, 3+7 *nTrk, 3+7 *nTrk, ROOT::Math::MatRepSym< double, 3+7 *nTrk > > &fitCov)
 

Private Attributes

VirtualKinematicParticleFactory pFactory
 
KinematicVertexFactory vFactory
 

Detailed Description

Class constructing te final output tree for the constrained vertex fitter. To be used by corresponding fitter only. Tree builders are scheduled for generalization: They should be inherited from the single generic class in the next version of the library.

Definition at line 15 of file ConstrainedTreeBuilderT.h.

Constructor & Destructor Documentation

◆ ConstrainedTreeBuilderT()

ConstrainedTreeBuilderT::ConstrainedTreeBuilderT ( )
inline

Definition at line 17 of file ConstrainedTreeBuilderT.h.

17 {}

◆ ~ConstrainedTreeBuilderT()

ConstrainedTreeBuilderT::~ConstrainedTreeBuilderT ( )
inline

Definition at line 19 of file ConstrainedTreeBuilderT.h.

19 {}

Member Function Documentation

◆ buildRealTree()

RefCountedKinematicTree ConstrainedTreeBuilderT::buildRealTree ( const RefCountedKinematicParticle  virtualParticle,
const RefCountedKinematicVertex  vtx,
const std::vector< RefCountedKinematicParticle > &  particles 
) const
private

Definition at line 5 of file ConstrainedTreeBuilderT.cc.

References ecalTrigSettings_cff::particles, KinematicVertexFactory::vertex(), vFactory, and L1BJetProducer_cff::vtx.

Referenced by buildTree().

8  {
9  //making a resulting tree:
11 
12  //fake production vertex:
14  resTree->addParticle(fVertex, vtx, virtualParticle);
15 
16  //adding final state
17  for (std::vector<RefCountedKinematicParticle>::const_iterator il = particles.begin(); il != particles.end(); il++) {
18  if ((*il)->previousParticle()->correspondingTree() != nullptr) {
19  KinematicTree* tree = (*il)->previousParticle()->correspondingTree();
20  tree->movePointerToTheTop();
21  tree->replaceCurrentParticle(*il);
22  RefCountedKinematicVertex cdVertex = resTree->currentDecayVertex();
23  resTree->addTree(cdVertex, tree);
24  } else {
26  resTree->addParticle(vtx, ffVertex, *il);
27  }
28  }
29  return resTree;
30 }
KinematicVertexFactory vFactory
static RefCountedKinematicVertex vertex(const VertexState &state, float totalChiSq, float degreesOfFr)
Definition: tree.py:1

◆ buildTree()

template<int nTrk>
RefCountedKinematicTree ConstrainedTreeBuilderT::buildTree ( const std::vector< RefCountedKinematicParticle > &  initialParticles,
const std::vector< KinematicState > &  finalStates,
const RefCountedKinematicVertex  vtx,
const ROOT::Math::SMatrix< double, 3+7 *nTrk, 3+7 *nTrk, ROOT::Math::MatRepSym< double, 3+7 *nTrk > > &  fCov 
) const

Method constructing tree out of set of refitted states, vertex, and full covariance matrix.

Definition at line 57 of file ConstrainedTreeBuilderT.h.

References a, buildRealTree(), ALCARECOTkAlJpsiMuMu_cff::charge, nano_mu_local_reco_cff::chi2, pseudoTop_cfi::finalStates, mps_fire::i, LogDebug, HLT_2023v12_cff::magneticField, dqmiodumpmetadata::n, hcaldqm::flag::nState, AlCaHLTBitMon_ParallelJobs::p, VirtualKinematicParticleFactory::particle(), pFactory, mathSSE::sqrt(), bphysicsOniaDQM_cfi::vertex, L1BJetProducer_cff::vtx, and z.

62  {
63  if (!vertex->vertexIsValid()) {
64  LogDebug("RecoVertex/ConstrainedTreeBuilder") << "Vertex is invalid\n";
66  }
68  vtx(0) = vertex->position().x();
69  vtx(1) = vertex->position().y();
70  vtx(2) = vertex->position().z();
71  AlgebraicMatrix33 vertexCov = fullCov.template Sub<ROOT::Math::SMatrix<double, 3> >(0, 0);
72 
73  // cout << fullCov<<endl;
74  // cout << "RecoVertex/ConstrainedTreeBuilder"<<vtx<<endl;
75 
76  double ent = 0.;
77  int charge = 0;
78  AlgebraicVector7 virtualPartPar;
79  virtualPartPar(0) = vertex->position().x();
80  virtualPartPar(1) = vertex->position().y();
81  virtualPartPar(2) = vertex->position().z();
82 
83  //making refitted particles out of refitted states.
84  //none of the operations above violates the order of particles
85 
86  ROOT::Math::SMatrix<double, 7, 7, ROOT::Math::MatRepStd<double, 7, 7> > aMatrix;
87  ROOT::Math::SMatrix<double, 7, 7, ROOT::Math::MatRepStd<double, 7, 7> > aMatrixT;
88  aMatrix(3, 3) = aMatrixT(3, 3) = 1;
89  aMatrix(4, 4) = aMatrixT(4, 4) = 1;
90  aMatrix(5, 5) = aMatrixT(5, 5) = 1;
91  aMatrix(6, 6) = aMatrixT(6, 6) = 1;
92  ROOT::Math::SMatrix<double, 7, 3, ROOT::Math::MatRepStd<double, 7, 3> > bMatrix;
93  ROOT::Math::SMatrix<double, 3, 7, ROOT::Math::MatRepStd<double, 3, 7> > bMatrixT;
94  bMatrix(0, 0) = bMatrixT(0, 0) = 1;
95  bMatrix(1, 1) = bMatrixT(1, 1) = 1;
96  bMatrix(2, 2) = bMatrixT(2, 2) = 1;
97  AlgebraicSymMatrix77 trackParCov;
98  ROOT::Math::SMatrix<double, 3, 7, ROOT::Math::MatRepStd<double, 3, 7> > vtxTrackCov;
99  AlgebraicSymMatrix77 nCovariance;
100  // AlgebraicSymMatrix77 tmp;
101 
102  std::vector<RefCountedKinematicParticle>::const_iterator i = initialParticles.begin();
103  std::vector<KinematicState>::const_iterator iStates = finalStates.begin();
104  std::vector<RefCountedKinematicParticle> rParticles;
105  int n = 0;
106  // assert(initialParticles.size()==nTrk);
107  for (; i != initialParticles.end() && iStates != finalStates.end(); ++i, ++iStates) {
108  AlgebraicVector7 p = iStates->kinematicParameters().vector();
109  double a = -iStates->particleCharge() * iStates->magneticField()->inInverseGeV(iStates->globalPosition()).z();
110 
111  aMatrix(4, 0) = aMatrixT(0, 4) = -a;
112  aMatrix(3, 1) = aMatrixT(1, 3) = a;
113  bMatrix(4, 0) = bMatrixT(0, 4) = a;
114  bMatrix(3, 1) = bMatrixT(1, 3) = -a;
115 
116  AlgebraicVector7 par = aMatrix * p + bMatrix * vtx;
117 
118  trackParCov = fullCov.template Sub<AlgebraicSymMatrix77>(3 + n * 7, 3 + n * 7);
119  vtxTrackCov = fullCov.template Sub<ROOT::Math::SMatrix<double, 3, 7> >(0, 3 + n * 7);
120  ROOT::Math::AssignSym::Evaluate(nCovariance,
121  aMatrix * trackParCov * aMatrixT +
122  aMatrix * ROOT::Math::Transpose(vtxTrackCov) * bMatrixT +
123  bMatrix * vtxTrackCov * aMatrixT + bMatrix * vertexCov * bMatrixT);
124  /*
125  ROOT::Math::AssignSym::Evaluate(tmp, aMatrix * ROOT::Math::Transpose(vtxTrackCov) * bMatrixT);
126  nCovariance+=tmp;
127  ROOT::Math::AssignSym::Evaluate(tmp, bMatrix * vtxTrackCov * aMatrixT);
128  nCovariance+=tmp;
129  ROOT::Math::AssignSym::Evaluate(tmp, bMatrix * vertexCov * bMatrixT);
130  nCovariance+=tmp;
131  */
132 
133  KinematicState stateAtVertex(KinematicParameters(par),
134  KinematicParametersError(nCovariance),
135  iStates->particleCharge(),
136  iStates->magneticField());
137  rParticles.push_back((*i)->refittedParticle(stateAtVertex, vertex->chiSquared(), vertex->degreesOfFreedom()));
138 
139  virtualPartPar(3) += par(3);
140  virtualPartPar(4) += par(4);
141  virtualPartPar(5) += par(5);
142  ent += sqrt(par(6) * par(6) + par(3) * par(3) + par(4) * par(4) + par(5) * par(5));
143  charge += iStates->particleCharge();
144 
145  ++n;
146  }
147 
148  //total reconstructed mass
149  double differ = ent * ent - (virtualPartPar(3) * virtualPartPar(3) + virtualPartPar(5) * virtualPartPar(5) +
150  virtualPartPar(4) * virtualPartPar(4));
151  if (differ > 0.) {
152  virtualPartPar(6) = sqrt(differ);
153  } else {
154  LogDebug("ConstrainedTreeBuilder") << "Fit failed: Current precision does not allow to calculate the mass\n";
156  }
157 
158  // covariance matrix:
159 
160  AlgebraicSymMatrix77 cov = this->covarianceMatrix<nTrk>(rParticles, virtualPartPar, fullCov);
161 
163  KinematicParameters(virtualPartPar), KinematicParametersError(cov), charge, initialParticles[0]->magneticField());
164 
165  //newborn kinematic particle
166  float chi2 = vertex->chiSquared();
167  float ndf = vertex->degreesOfFreedom();
168  KinematicParticle* zp = nullptr;
169  RefCountedKinematicParticle virtualParticle = pFactory.particle(nState, chi2, ndf, zp);
170 
171  return buildRealTree(virtualParticle, vertex, rParticles);
172 }
ROOT::Math::SMatrix< double, 3, 3, ROOT::Math::MatRepStd< double, 3, 3 > > AlgebraicMatrix33
VirtualKinematicParticleFactory pFactory
ROOT::Math::SVector< double, 7 > AlgebraicVector7
Definition: Matrices.h:8
T sqrt(T t)
Definition: SSEVec.h:19
RefCountedKinematicTree buildRealTree(const RefCountedKinematicParticle virtualParticle, const RefCountedKinematicVertex vtx, const std::vector< RefCountedKinematicParticle > &particles) const
ROOT::Math::SMatrix< double, 7, 7, ROOT::Math::MatRepSym< double, 7 > > AlgebraicSymMatrix77
Definition: Matrices.h:9
RefCountedKinematicParticle particle(const KinematicState &kineState, float &chiSquared, float &degreesOfFr, ReferenceCountingPointer< KinematicParticle > previousParticle, KinematicConstraint *lastConstraint=nullptr) const
double a
Definition: hdecay.h:121
ROOT::Math::SVector< double, 3 > AlgebraicVector3
#define LogDebug(id)

◆ covarianceMatrix()

template<int nTrk>
AlgebraicSymMatrix77 ConstrainedTreeBuilderT::covarianceMatrix ( const std::vector< RefCountedKinematicParticle > &  rPart,
const AlgebraicVector7 newPar,
const ROOT::Math::SMatrix< double, 3+7 *nTrk, 3+7 *nTrk, ROOT::Math::MatRepSym< double, 3+7 *nTrk > > &  fitCov 
)
staticprivate

Metod to reconstruct the full covariance matrix of the resulting particle.

Definition at line 175 of file ConstrainedTreeBuilderT.h.

References mps_fire::i, dqmiolumiharvest::j, HLT_2023v12_cff::magneticField, runTheMatrix::ret, findQualityFiles::size, mathSSE::sqrt(), and HcalDetIdTransform::transform().

178  {
179  typedef ROOT::Math::SMatrix<double, 3 + 7 * nTrk, 3 + 7 * nTrk, ROOT::Math::MatRepSym<double, 3 + 7 * nTrk> > FitCov;
180  //constructing the full matrix using the simple fact
181  //that we have never broken the order of tracks
182  // during our fit.
183 
184  int size = nTrk;
185 
186  if (int(rPart.size()) != size)
187  throw "error in ConstrainedTreeBuilderT ";
188 
189  //global propagation to the vertex position
190  //Jacobian is done for all the parameters together
191  ROOT::Math::SMatrix<double, 3 + 7 * nTrk> jac;
192  jac(0, 0) = 1;
193  jac(1, 1) = 1;
194  jac(2, 2) = 1;
195  ROOT::Math::SMatrix<double, 3, 7> upper;
196  ROOT::Math::SMatrix<double, 7> diagonal;
197  for (int i_int = 0; i_int != size; ++i_int) {
198  RefCountedKinematicParticle const& i = rPart[i_int];
199  //vertex position related components of the matrix
200  double a_i = -(i)->currentState().particleCharge() *
201  (i)->magneticField()->inInverseGeV((i)->currentState().globalPosition()).z();
202  upper(0, 0) = 1;
203  upper(1, 1) = 1;
204  upper(2, 2) = 1;
205  upper(1, 3) = -a_i;
206  upper(0, 4) = a_i;
207  jac.Place_at(upper, 0, 3 + i_int * 7);
208 
209  diagonal(3, 3) = 1;
210  diagonal(4, 4) = 1;
211  diagonal(5, 5) = 1;
212  diagonal(6, 6) = 1;
213  diagonal(1, 3) = a_i;
214  diagonal(0, 4) = -a_i;
215  jac.Place_at(diagonal, 3 + i_int * 7, 3 + i_int * 7);
216  }
217 
218  // jacobian is constructed in such a way, that
219  // right operation for transformation will be
220  // fitCov.similarityT(jac)
221  // WARNING: normal similarity operation is
222  // not valid in this case
223  //now making reduced matrix:
224  // int vSize = rPart.size();
225  FitCov const& fit_cov_sym = fitCov;
226  /*
227  for(int i = 0; i<7*vSize+3; ++i)
228  {
229  for(int j = 0; j<7*vSize+3; ++j)
230  {if(i<=j) fit_cov_sym(i,j) = fitCov(i,j);}
231  }
232  */
233 
234  ROOT::Math::SMatrix<double, 4 * nTrk + 3, 4 * nTrk + 3, ROOT::Math::MatRepSym<double, 4 * nTrk + 3> > reduced;
235  FitCov transform = ROOT::Math::SimilarityT(jac, fit_cov_sym); // similarityT???
236 
237  //jacobian to add matrix components
238  ROOT::Math::SMatrix<double, 7, 4 * nTrk + 3> jac_t;
239  jac_t(0, 0) = 1.;
240  jac_t(1, 1) = 1.;
241  jac_t(2, 2) = 1.;
242 
243  double energy_global =
244  sqrt(newPar(3) * newPar(3) + newPar(4) * newPar(4) + newPar(5) * newPar(5) + newPar(6) * newPar(6));
245  for (int il_int = 0; il_int != size; ++il_int) {
246  RefCountedKinematicParticle const& rs = rPart[il_int];
247  //jacobian components:
248  int off1 = 3;
249  int off2 = il_int * 4 + 3;
250  jac_t(off1 + 0, off2 + 0) = 1.;
251  jac_t(off1 + 1, off2 + 1) = 1.;
252  jac_t(off1 + 2, off2 + 2) = 1.;
253 
254  //non-trival elements: mass correlations:
255  AlgebraicVector7 l_Par = (rs)->currentState().kinematicParameters().vector();
256  double energy_local = sqrt(l_Par(6) * l_Par(6) + l_Par(3) * l_Par(3) + l_Par(4) * l_Par(4) + l_Par(5) * l_Par(5));
257  jac_t(off1 + 3, off2 + 3) = energy_global * l_Par(6) / (newPar(6) * energy_local);
258  jac_t(off1 + 3, off2 + 0) = ((energy_global * l_Par(3) / energy_local) - newPar(3)) / newPar(6);
259  jac_t(off1 + 3, off2 + 1) = ((energy_global * l_Par(4) / energy_local) - newPar(4)) / newPar(6);
260  jac_t(off1 + 3, off2 + 2) = ((energy_global * l_Par(5) / energy_local) - newPar(5)) / newPar(6);
261  }
262 
263  for (int i = 0; i < 7; ++i)
264  for (int j = 0; j < 7; ++j)
265  reduced(i, j) = transform(i + 3, j + 3);
266 
267  for (int i = 1; i < size; i++) {
268  //non-trival elements: mass correlations:
269  //upper row and right column
270 
271  int off1 = 0;
272  int off2 = 3 + 4 * i;
273  for (int l1 = 0; l1 < 3; ++l1)
274  for (int l2 = 0; l2 < 4; ++l2)
275  reduced(off1 + l1, off2 + l2) = transform(3 + l1, 6 + 7 * i + l2);
276 
277  //diagonal elements
278  off1 = off2 = 3 + 4 * i;
279  for (int l1 = 0; l1 < 4; ++l1)
280  for (int l2 = 0; l2 < 4; ++l2)
281  reduced(off1 + l1, off2 + l2) = transform(6 + 7 * i + l1, 6 + 7 * i + l2);
282 
283  //off diagonal elements
284  for (int j = 1; j < size; j++) {
285  off1 = 3 + 4 * (i - 1);
286  off2 = 3 + 4 * j;
287  for (int l1 = 0; l1 < 4; ++l1)
288  for (int l2 = 0; l2 < 4; ++l2)
289  reduced(off1 + l1, off2 + l2) = transform(6 + 7 * (i - 1) + l1, 6 + 7 * j + l2);
290  }
291  }
292 
294  ROOT::Math::AssignSym::Evaluate(ret, jac_t * reduced * ROOT::Math::Transpose(jac_t));
295  return ret;
296 }
size
Write out results.
ret
prodAgent to be discontinued
ROOT::Math::SVector< double, 7 > AlgebraicVector7
Definition: Matrices.h:8
T sqrt(T t)
Definition: SSEVec.h:19
ROOT::Math::SMatrix< double, 7, 7, ROOT::Math::MatRepSym< double, 7 > > AlgebraicSymMatrix77
Definition: Matrices.h:9
unsigned transform(const HcalDetId &id, unsigned transformCode)

Member Data Documentation

◆ pFactory

VirtualKinematicParticleFactory ConstrainedTreeBuilderT::pFactory
private

Definition at line 49 of file ConstrainedTreeBuilderT.h.

Referenced by buildTree().

◆ vFactory

KinematicVertexFactory ConstrainedTreeBuilderT::vFactory
private

Definition at line 50 of file ConstrainedTreeBuilderT.h.

Referenced by buildRealTree().