CMS 3D CMS Logo

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

#include <FinalTreeBuilder.h>

Public Member Functions

RefCountedKinematicTree buildTree (const CachingVertex< 6 > &vtx, const std::vector< RefCountedKinematicParticle > &input) const
 
 FinalTreeBuilder ()
 
 ~FinalTreeBuilder ()
 

Private Types

typedef ReferenceCountingPointer< LinearizedTrackState< 6 > > RefCountedLinearizedTrackState
 
typedef ReferenceCountingPointer< RefittedTrackState< 6 > > RefCountedRefittedTrackState
 
typedef ReferenceCountingPointer< VertexTrack< 6 > > RefCountedVertexTrack
 

Private Member Functions

AlgebraicMatrix momentumPart (const CachingVertex< 6 > &vtx, const AlgebraicVector7 &par) const
 

Private Attributes

KinematicVertexFactorykvFactory
 
VirtualKinematicParticleFactorypFactory
 

Detailed Description

Class building a resulting output tree out of the information provided by KinematicParticleVertexFitter.

Definition at line 16 of file FinalTreeBuilder.h.

Member Typedef Documentation

◆ RefCountedLinearizedTrackState

Definition at line 27 of file FinalTreeBuilder.h.

◆ RefCountedRefittedTrackState

Definition at line 28 of file FinalTreeBuilder.h.

◆ RefCountedVertexTrack

Definition at line 26 of file FinalTreeBuilder.h.

Constructor & Destructor Documentation

◆ FinalTreeBuilder()

FinalTreeBuilder::FinalTreeBuilder ( )

◆ ~FinalTreeBuilder()

FinalTreeBuilder::~FinalTreeBuilder ( )

Definition at line 14 of file FinalTreeBuilder.cc.

References kvFactory, and pFactory.

14  {
15  delete kvFactory;
16  delete pFactory;
17 }
VirtualKinematicParticleFactory * pFactory
KinematicVertexFactory * kvFactory

Member Function Documentation

◆ buildTree()

RefCountedKinematicTree FinalTreeBuilder::buildTree ( const CachingVertex< 6 > &  vtx,
const std::vector< RefCountedKinematicParticle > &  input 
) const

Definition at line 19 of file FinalTreeBuilder.cc.

References mps_fire::i, input, dqmiolumiharvest::j, KinematicRefittedTrackState::kinematicMomentumVector(), KinematicRefittedTrackState::kinematicParameters(), KinematicRefittedTrackState::kinematicParametersCovariance(), kvFactory, LogDebug, momentumPart(), nPart(), hcaldqm::flag::nState, VirtualKinematicParticleFactory::particle(), pFactory, mathSSE::sqrt(), KinematicVertexFactory::vertex(), and L1BJetProducer_cff::vtx.

Referenced by KinematicParticleVertexFitter::fit().

20  {
21  //creating resulting kinematic particle
22  AlgebraicVector7 par;
23  AlgebraicMatrix cov(7, 7, 0);
24  par(0) = vtx.position().x();
25  par(1) = vtx.position().y();
26  par(2) = vtx.position().z();
27  double en = 0.;
28  int ch = 0;
29 
30  //new particle momentum calculation and refitted kinematic states
31  std::vector<KinematicRefittedTrackState*> rStates;
32  std::vector<RefCountedVertexTrack> refTracks = vtx.tracks();
33  for (std::vector<RefCountedVertexTrack>::const_iterator i = refTracks.begin(); i != refTracks.end(); ++i) {
34  KinematicRefittedTrackState* rs = dynamic_cast<KinematicRefittedTrackState*>(&(*((*i)->refittedState())));
36  par(3) += f_mom(0);
37  par(4) += f_mom(1);
38  par(5) += f_mom(2);
39  en += sqrt(f_mom(1) * f_mom(1) + f_mom(2) * f_mom(2) + f_mom(3) * f_mom(3) + f_mom(0) * f_mom(0));
40  ch += (*i)->linearizedTrack()->charge();
41  rStates.push_back(rs);
42  }
43 
44  //math precision check (numerical stability)
45  double differ = en * en - (par(3) * par(3) + par(4) * par(4) + par(5) * par(5));
46  if (differ > 0.) {
47  par(6) = sqrt(differ);
48  } else {
49  LogDebug("FinalTreeBuilder") << "Fit failed: Current precision does not allow to calculate the mass\n";
51  }
52 
53  // covariance matrix calculation: momentum-momentum components part (4x4)
54  // and position-momentum components part:
55  AlgebraicMatrix m_all = momentumPart(vtx, par);
56 
57  //position-position components part (3x3)
58  // AlgebraicMatrix x_X = vtx.error().matrix();
59 
60  //making new matrix itself
61  cov.sub(1, 1, m_all);
62 
63  //covariance sym matrix
65  for (int i = 1; i < 8; i++) {
66  for (int j = 1; j < 8; j++) {
67  if (i <= j)
68  sCov(i - 1, j - 1) = cov(i, j);
69  }
70  }
71 
72  //valid decay vertex for our resulting particle
74 
75  //invalid production vertex for our resulting particle
77 
78  //new born kinematic particle
79  KinematicParameters kPar(par);
80  KinematicParametersError kEr(sCov);
81  const MagneticField* field = input.front()->magneticField();
82  KinematicState nState(kPar, kEr, ch, field);
83 
84  //invalid previous particle and empty constraint:
85  KinematicParticle* zp = nullptr;
87 
88  float vChi = vtx.totalChiSquared();
89  float vNdf = vtx.degreesOfFreedom();
91 
92  //adding top particle to the tree
94  resTree->addParticle(pVrt, dVrt, nPart);
95 
96  //making refitted kinematic particles and putting them to the tree
97  std::vector<RefCountedKinematicParticle> rrP;
98 
99  std::vector<RefCountedKinematicParticle>::const_iterator j;
100  std::vector<RefCountedVertexTrack>::const_iterator i;
101  for (j = input.begin(), i = refTracks.begin(); j != input.end() && i != refTracks.end(); ++j, ++i) {
102  RefCountedLinearizedTrackState lT = (*i)->linearizedTrack();
103  KinematicRefittedTrackState* rS = dynamic_cast<KinematicRefittedTrackState*>(&(*((*i)->refittedState())));
104 
105  // RefCountedRefittedTrackState rS = (*i)->refittedState();
107  KinematicParameters lkPar(lPar);
109  KinematicParametersError lkCov(lCov);
110  TrackCharge lch = lT->charge();
111  KinematicState nState(lkPar, lkCov, lch, field);
112  RefCountedKinematicParticle nPart = (*j)->refittedParticle(nState, vChi, vNdf);
113  rrP.push_back(nPart);
114  if ((*j)->correspondingTree() != nullptr) {
115  //here are the particles having trees after them
116  KinematicTree* tree = (*j)->correspondingTree();
117  tree->movePointerToTheTop();
118  tree->replaceCurrentParticle(nPart);
119  RefCountedKinematicVertex cdVertex = resTree->currentDecayVertex();
120  resTree->addTree(cdVertex, tree);
121  } else {
122  //here are just particles fitted to this tree
124  resTree->addParticle(dVrt, nV, nPart);
125  }
126  }
127  return resTree;
128 }
ReferenceCountingPointer< LinearizedTrackState< 6 > > RefCountedLinearizedTrackState
ROOT::Math::SVector< double, 7 > AlgebraicVector7
Definition: Matrices.h:8
VirtualKinematicParticleFactory * pFactory
static std::string const input
Definition: EdmProvDump.cc:50
AlgebraicVector4 kinematicMomentumVector() const
int TrackCharge
Definition: TrackCharge.h:4
CLHEP::HepMatrix AlgebraicMatrix
T sqrt(T t)
Definition: SSEVec.h:19
ROOT::Math::SVector< double, 4 > AlgebraicVector4
KinematicVertexFactory * kvFactory
TString nPart(Int_t part, TString string, TString delimit=";", Bool_t removerest=true)
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
AlgebraicVector7 kinematicParameters() const
AlgebraicSymMatrix77 kinematicParametersCovariance() const
static RefCountedKinematicVertex vertex(const VertexState &state, float totalChiSq, float degreesOfFr)
Definition: tree.py:1
AlgebraicMatrix momentumPart(const CachingVertex< 6 > &vtx, const AlgebraicVector7 &par) const
#define LogDebug(id)

◆ momentumPart()

AlgebraicMatrix FinalTreeBuilder::momentumPart ( const CachingVertex< 6 > &  vtx,
const AlgebraicVector7 par 
) const
private

Definition at line 132 of file FinalTreeBuilder.cc.

References a, funct::cos(), Exception, EgHLTOffHistBins_cfi::mass, phi, rho, funct::sin(), findQualityFiles::size, mathSSE::sqrt(), funct::tan(), theta(), L1BJetProducer_cff::vtx, and z.

Referenced by buildTree().

132  {
133  std::vector<RefCountedVertexTrack> refTracks = vtx.tracks();
134  int size = refTracks.size();
135  AlgebraicMatrix cov(7 + 4 * (size - 1), 7 + 4 * (size - 1));
136  AlgebraicMatrix jac(7, 7 + 4 * (size - 1));
137  jac(1, 1) = 1.;
138  jac(2, 2) = 1.;
139  jac(3, 3) = 1.;
140 
141  double energy_total = sqrt(par(3) * par(3) + par(6) * par(6) + par(5) * par(5) + par(4) * par(4));
142 
143  std::vector<RefCountedVertexTrack>::const_iterator rt_i;
144  int i_int = 0;
145  for (rt_i = refTracks.begin(); rt_i != refTracks.end(); rt_i++) {
146  double a;
147  AlgebraicVector6 param = (**rt_i).refittedState()->parameters(); // rho, theta, phi,tr_im, z_im, mass
148  double rho = param[0];
149  double theta = param[1];
150  double phi = param[2];
151  double mass = param[5];
152 
153  if ((**rt_i).linearizedTrack()->charge() != 0) {
154  a = -(**rt_i).refittedState()->freeTrajectoryState().parameters().magneticFieldInInverseGeV(vtx.position()).z() *
155  (**rt_i).refittedState()->freeTrajectoryState().parameters().charge();
156  if (a == 0.)
157  throw cms::Exception("FinalTreeBuilder", "Field is 0");
158  } else {
159  a = 1;
160  }
161 
162  AlgebraicMatrix jc_el(4, 4, 0);
163  jc_el(1, 1) = -a * cos(phi) / (rho * rho); //dpx/d rho
164  jc_el(2, 1) = -a * sin(phi) / (rho * rho); //dpy/d rho
165  jc_el(3, 1) = -a / (rho * rho * tan(theta)); //dpz/d rho
166 
167  jc_el(3, 2) = -a / (rho * sin(theta) * sin(theta)); //dpz/d theta
168 
169  jc_el(1, 3) = -a * sin(phi) / rho; //dpx/d phi
170  jc_el(2, 3) = a * cos(phi) / rho; //dpy/d
171 
172  //non-trival elements: mass correlations:
173  double energy_local = sqrt(a * a / (rho * rho) * (1 + 1 / (tan(theta) * tan(theta))) + mass * mass);
174 
175  jc_el(4, 4) = energy_total * mass / (par(6) * energy_local); // dm/dm
176 
177  jc_el(4, 1) = (-(energy_total / energy_local * a * a / (rho * rho * rho * sin(theta) * sin(theta))) +
178  par(3) * a / (rho * rho) * cos(phi) + par(4) * a / (rho * rho) * sin(phi) +
179  par(5) * a / (rho * rho * tan(theta))) /
180  par(6); //dm / drho
181 
182  jc_el(4, 2) = (-(energy_total / energy_local * a * a / (rho * rho * sin(theta) * sin(theta) * tan(theta))) +
183  par(5) * a / (rho * sin(theta) * sin(theta))) /
184  par(6); //dm d theta
185 
186  jc_el(4, 3) = (par(3) * sin(phi) - par(4) * cos(phi)) * a / (rho * par(6)); //dm/dphi
187 
188  jac.sub(4, i_int * 4 + 4, jc_el);
189 
190  //top left corner elements
191  if (i_int == 0) {
192  cov.sub(1, 1, asHepMatrix<7>((**rt_i).fullCovariance()));
193  } else {
194  //4-momentum corellatons: diagonal elements of the matrix
195  AlgebraicMatrix fullCovMatrix(asHepMatrix<7>((**rt_i).fullCovariance()));
196  AlgebraicMatrix m_m_cov = fullCovMatrix.sub(4, 7, 4, 7);
197  AlgebraicMatrix x_p_cov = fullCovMatrix.sub(1, 3, 4, 7);
198  AlgebraicMatrix p_x_cov = fullCovMatrix.sub(4, 7, 1, 3);
199 
200  // cout << "Full covariance: \n"<< (**rt_i).fullCovariance()<<endl;
201  // cout << "Full m_m_cov: "<< m_m_cov<<endl;
202  // cout<<"p_x_cov"<< p_x_cov<<endl;
203  // cout<<"x_p_cov"<< x_p_cov<<endl;
204 
205  //putting everything to the joint covariance matrix:
206  //diagonal momentum-momentum elements:
207  cov.sub(i_int * 4 + 4, i_int * 4 + 4, m_m_cov);
208 
209  //position momentum elements:
210  cov.sub(1, i_int * 4 + 4, x_p_cov);
211  cov.sub(i_int * 4 + 4, 1, p_x_cov);
212 
213  //off diagonal elements: corellations
214  // track momentum - track momentum
215  }
216  int j_int = 0;
217  for (std::vector<RefCountedVertexTrack>::const_iterator rt_j = refTracks.begin(); rt_j != refTracks.end(); rt_j++) {
218  if (i_int < j_int) {
219  AlgebraicMatrix i_k_cov_m = asHepMatrix<4, 4>(vtx.tkToTkCovariance((*rt_i), (*rt_j)));
220  // cout<<"i_k_cov_m"<<i_k_cov_m <<endl;
221  cov.sub(i_int * 4 + 4, j_int * 4 + 4, i_k_cov_m);
222  cov.sub(j_int * 4 + 4, i_int * 4 + 4, i_k_cov_m.T());
223  }
224  j_int++;
225  }
226  i_int++;
227  }
228  // cout<<"jac"<<jac<<endl;
229  // cout<<"cov"<<cov<<endl;
230  // cout << "final result new"<<jac*cov*jac.T()<<endl;
231 
232  return jac * cov * jac.T();
233 }
size
Write out results.
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
CLHEP::HepMatrix AlgebraicMatrix
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
double a
Definition: hdecay.h:121
ROOT::Math::SVector< double, 6 > AlgebraicVector6
Geom::Theta< T > theta() const

Member Data Documentation

◆ kvFactory

KinematicVertexFactory* FinalTreeBuilder::kvFactory
private

Definition at line 33 of file FinalTreeBuilder.h.

Referenced by buildTree(), FinalTreeBuilder(), and ~FinalTreeBuilder().

◆ pFactory

VirtualKinematicParticleFactory* FinalTreeBuilder::pFactory
private

Definition at line 34 of file FinalTreeBuilder.h.

Referenced by buildTree(), FinalTreeBuilder(), and ~FinalTreeBuilder().