CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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

Definition at line 29 of file FinalTreeBuilder.h.

Definition at line 30 of file FinalTreeBuilder.h.

Definition at line 28 of file FinalTreeBuilder.h.

Constructor & Destructor Documentation

FinalTreeBuilder::FinalTreeBuilder ( )
FinalTreeBuilder::~FinalTreeBuilder ( )

Definition at line 15 of file FinalTreeBuilder.cc.

References kvFactory, and pFactory.

16 {
17  delete kvFactory;
18  delete pFactory;
19 }
VirtualKinematicParticleFactory * pFactory
KinematicVertexFactory * kvFactory

Member Function Documentation

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

Definition at line 21 of file FinalTreeBuilder.cc.

References CachingVertex< N >::degreesOfFreedom(), i, j, KinematicRefittedTrackState::kinematicMomentumVector(), KinematicRefittedTrackState::kinematicParameters(), KinematicRefittedTrackState::kinematicParametersCovariance(), kvFactory, LogDebug, momentumPart(), KinematicTree::movePointerToTheTop(), nPart(), VirtualKinematicParticleFactory::particle(), pFactory, CachingVertex< N >::position(), KinematicTree::replaceCurrentParticle(), mathSSE::sqrt(), CachingVertex< N >::totalChiSquared(), CachingVertex< N >::tracks(), MainPageGenerator::tree, KinematicVertexFactory::vertex(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by KinematicParticleVertexFitter::fit().

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

Definition at line 143 of file FinalTreeBuilder.cc.

References a, funct::cos(), Exception, phi, CachingVertex< N >::position(), rho, funct::sin(), findQualityFiles::size, mathSSE::sqrt(), funct::tan(), theta(), CachingVertex< N >::tkToTkCovariance(), CachingVertex< N >::tracks(), and z.

Referenced by buildTree().

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

Member Data Documentation

KinematicVertexFactory* FinalTreeBuilder::kvFactory
private

Definition at line 36 of file FinalTreeBuilder.h.

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

VirtualKinematicParticleFactory* FinalTreeBuilder::pFactory
private

Definition at line 37 of file FinalTreeBuilder.h.

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