#include <RecoVertex/KinematicFit/interface/FinalTreeBuilder.h>
Public Member Functions | |
RefCountedKinematicTree | buildTree (const CachingVertex< 6 > &vtx, 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 (vector< KinematicRefittedTrackState * > rStates, const CachingVertex< 6 > &vtx, const AlgebraicVector7 &par) const |
Private Attributes | |
KinematicVertexFactory * | kvFactory |
VirtualKinematicParticleFactory * | pFactory |
Definition at line 16 of file FinalTreeBuilder.h.
typedef ReferenceCountingPointer<LinearizedTrackState<6> > FinalTreeBuilder::RefCountedLinearizedTrackState [private] |
Definition at line 29 of file FinalTreeBuilder.h.
typedef ReferenceCountingPointer<RefittedTrackState<6> > FinalTreeBuilder::RefCountedRefittedTrackState [private] |
Definition at line 30 of file FinalTreeBuilder.h.
typedef ReferenceCountingPointer<VertexTrack<6> > FinalTreeBuilder::RefCountedVertexTrack [private] |
Definition at line 28 of file FinalTreeBuilder.h.
FinalTreeBuilder::FinalTreeBuilder | ( | ) |
Definition at line 8 of file FinalTreeBuilder.cc.
References kvFactory, and pFactory.
00009 { 00010 kvFactory = new KinematicVertexFactory(); 00011 KinematicStatePropagator * ksp = 0; 00012 pFactory = new VirtualKinematicParticleFactory(ksp); 00013 }
FinalTreeBuilder::~FinalTreeBuilder | ( | ) |
RefCountedKinematicTree FinalTreeBuilder::buildTree | ( | const CachingVertex< 6 > & | vtx, | |
vector< RefCountedKinematicParticle > | input | |||
) | const |
Definition at line 21 of file FinalTreeBuilder.cc.
References GenMuonPlsPt100GeV_cfg::cout, CachingVertex< N >::degreesOfFreedom(), lat::endl(), CachingVertex< N >::error(), i, j, KinematicRefittedTrackState::kinematicMomentumVector(), KinematicRefittedTrackState::kinematicParameters(), KinematicRefittedTrackState::kinematicParametersCovariance(), kvFactory, GlobalErrorBase< T, ErrorWeightType >::matrix(), momentumPart(), KinematicTree::movePointerToTheTop(), VirtualKinematicParticleFactory::particle(), pFactory, CachingVertex< N >::position(), KinematicTree::replaceCurrentParticle(), funct::sqrt(), CachingVertex< N >::totalChiSquared(), CachingVertex< N >::tracks(), tree, KinematicVertexFactory::vertex(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().
Referenced by KinematicParticleVertexFitter::fit().
00023 { 00024 //creating resulting kinematic particle 00025 AlgebraicVector7 par; 00026 AlgebraicMatrix cov(7,7,0); 00027 par(0) = vtx.position().x(); 00028 par(1) = vtx.position().y(); 00029 par(2) = vtx.position().z(); 00030 double en = 0.; 00031 int ch = 0; 00032 00033 //new particle momentum calculation and refitted kinematic states 00034 vector<KinematicRefittedTrackState *> rStates; 00035 vector<RefCountedVertexTrack> refTracks = vtx.tracks(); 00036 for(vector<RefCountedVertexTrack>::const_iterator i = refTracks.begin();i !=refTracks.end();++i) 00037 { 00038 KinematicRefittedTrackState * rs = dynamic_cast<KinematicRefittedTrackState *>(&(*((*i)->refittedState()))); 00039 AlgebraicVector4 f_mom = rs->kinematicMomentumVector(); 00040 par(3) += f_mom(0); 00041 par(4) += f_mom(1); 00042 par(5) += f_mom(2); 00043 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)); 00044 ch += (*i)->linearizedTrack()->charge(); 00045 rStates.push_back(rs); 00046 } 00047 00048 //math precision check (numerical stability) 00049 double differ = en*en - (par(3)*par(3)+par(4)*par(4)+par(5)*par(5)); 00050 if(differ>0.) 00051 { 00052 par(6) = sqrt(differ); 00053 }else{ 00054 cout<<"warning! current precision does not allow to calculate the mass"<<endl; 00055 throw VertexException("warning! current precision does not allow to calculate the mass"); 00056 00057 par(6) = 0.; 00058 } 00059 00060 // covariance matrix calculation: momentum-momentum components part (4x4) 00061 // and position-momentum components part: 00062 AlgebraicMatrix m_all = momentumPart(rStates,vtx,par); 00063 00064 //position-position components part (3x3) 00065 AlgebraicMatrix x_X = vtx.error().matrix(); 00066 00067 //making new matrix itself 00068 cov.sub(1,1,m_all); 00069 00070 //covariance sym matrix 00071 AlgebraicSymMatrix77 sCov; 00072 for(int i = 1; i<8; i++) 00073 { 00074 for(int j = 1; j<8; j++) 00075 { 00076 if(i<=j) sCov(i-1,j-1) = cov(i,j); 00077 } 00078 } 00079 00080 //valid decay vertex for our resulting particle 00081 RefCountedKinematicVertex dVrt = kvFactory->vertex(vtx); 00082 00083 //invalid production vertex for our resulting particle 00084 RefCountedKinematicVertex pVrt = kvFactory->vertex(); 00085 00086 //new born kinematic particle 00087 KinematicParameters kPar(par); 00088 KinematicParametersError kEr(sCov); 00089 const MagneticField* field=input.front()->magneticField(); 00090 KinematicState nState(kPar, kEr, ch, field); 00091 00092 //invalid previous particle and empty constraint: 00093 KinematicParticle * zp = 0; 00094 RefCountedKinematicParticle pPart = ReferenceCountingPointer<KinematicParticle>(zp); 00095 00096 float vChi = vtx.totalChiSquared(); 00097 float vNdf = vtx.degreesOfFreedom(); 00098 RefCountedKinematicParticle nPart = pFactory->particle(nState, vChi, vNdf, pPart); 00099 00100 //adding top particle to the tree 00101 RefCountedKinematicTree resTree = ReferenceCountingPointer<KinematicTree>(new KinematicTree()); 00102 resTree->addParticle(pVrt,dVrt,nPart); 00103 00104 //making refitted kinematic particles and putting them to the tree 00105 vector<RefCountedKinematicParticle> rrP; 00106 00107 vector<RefCountedKinematicParticle>::const_iterator j; 00108 vector<RefCountedVertexTrack>::const_iterator i; 00109 for(j=input.begin(), i=refTracks.begin(); j !=input.end(), i !=refTracks.end();++j, ++i) 00110 { 00111 RefCountedLinearizedTrackState lT = (*i)->linearizedTrack(); 00112 KinematicRefittedTrackState * rS= dynamic_cast<KinematicRefittedTrackState *>(&(*((*i)->refittedState()))); 00113 00114 // RefCountedRefittedTrackState rS = (*i)->refittedState(); 00115 AlgebraicVector7 lPar = rS->kinematicParameters(); 00116 KinematicParameters lkPar(lPar); 00117 AlgebraicSymMatrix77 lCov = rS->kinematicParametersCovariance(); 00118 KinematicParametersError lkCov(lCov); 00119 TrackCharge lch = lT->charge(); 00120 KinematicState nState(lkPar,lkCov,lch, field); 00121 RefCountedKinematicParticle nPart = (*j)->refittedParticle(nState,vChi,vNdf); 00122 rrP.push_back(nPart); 00123 if((*j)->correspondingTree() != 0) 00124 { 00125 00126 //here are the particles having trees after them 00127 KinematicTree * tree = (*j)->correspondingTree(); 00128 tree->movePointerToTheTop(); 00129 tree->replaceCurrentParticle(nPart); 00130 RefCountedKinematicVertex cdVertex = resTree->currentDecayVertex(); 00131 resTree->addTree(cdVertex, tree); 00132 }else{ 00133 00134 //here are just particles fitted to this tree 00135 RefCountedKinematicVertex nV = kvFactory->vertex(); 00136 resTree->addParticle(dVrt,nV,nPart); 00137 } 00138 } 00139 return resTree; 00140 }
AlgebraicMatrix FinalTreeBuilder::momentumPart | ( | vector< KinematicRefittedTrackState * > | rStates, | |
const CachingVertex< 6 > & | vtx, | |||
const AlgebraicVector7 & | par | |||
) | const [private] |
Definition at line 144 of file FinalTreeBuilder.cc.
References l1, l2, size, funct::sqrt(), CachingVertex< N >::tkToTkCovariance(), and CachingVertex< N >::tracks().
Referenced by buildTree().
00146 { 00147 vector<RefCountedVertexTrack> refTracks = vtx.tracks(); 00148 int size = rStates.size(); 00149 AlgebraicMatrix cov(7+4*(size-1),7+4*(size-1)); 00150 AlgebraicMatrix jac(7,7+4*(size-1)); 00151 jac(1,1) = 1.; 00152 jac(2,2) = 1.; 00153 jac(3,3) = 1.; 00154 vector<KinematicRefittedTrackState *>::const_iterator rs; 00155 vector<RefCountedVertexTrack>::const_iterator rt_i; 00156 int i_int = 0; 00157 for(rs = rStates.begin(), rt_i = refTracks.begin(); rs != rStates.end(), rt_i != refTracks.end(); rs++, rt_i++) 00158 { 00159 AlgebraicMatrix jc_el(4,4,0); 00160 jc_el(1,1) = 1.; 00161 jc_el(2,2) = 1.; 00162 jc_el(3,3) = 1.; 00163 00164 //non-trival elements: mass correlations: 00165 AlgebraicVector7 l_Par = (*rs)->kinematicParameters(); 00166 double energy_local = sqrt(l_Par(3)*l_Par(3) + l_Par(4)*l_Par(4) + l_Par(5)*l_Par(5) + l_Par(6)*l_Par(6)); 00167 00168 double energy_global = sqrt(par(3)*par(3)+par(6)*par(6) + par(5)*par(5)+par(4)*par(4)); 00169 00170 jc_el(4,4) = energy_global*l_Par(6)/(par(6)*energy_local); 00171 00172 jc_el(4,1) = ((energy_global*l_Par(3)/energy_local) - par(3))/par(6); 00173 jc_el(4,2) = ((energy_global*l_Par(4)/energy_local) - par(4))/par(6); 00174 jc_el(4,3) = ((energy_global*l_Par(5)/energy_local) - par(5))/par(6); 00175 00176 jac.sub(4,i_int*4+4,jc_el); 00177 00178 //top left corner elements 00179 if(i_int == 0) 00180 { 00181 cov.sub(1,1,asHepMatrix<7>((*rs)->kinematicParametersCovariance()) ); 00182 }else{ 00183 //4-momentum corellatons: diagonal elements of the matrix 00184 AlgebraicMatrix m_m_cov = asHepMatrix<7>((*rs)->kinematicParametersCovariance()).sub(4,7); 00185 00186 //position momentum and momentum position corellations 00187 AlgebraicMatrix xpcov = asHepMatrix<7>((*rs)->kinematicParametersCovariance()); 00188 AlgebraicMatrix p_x_cov(4,3); 00189 AlgebraicMatrix x_p_cov(3,4); 00190 00191 for(int l1 = 1; l1<5; ++l1) 00192 { 00193 for(int l2 = 1; l2<4; ++l2) 00194 {p_x_cov(l1,l2) = xpcov(3+l1,l2);} 00195 } 00196 00197 for(int l1 = 1; l1<4; ++l1) 00198 { 00199 for(int l2 = 1; l2<5; ++l2) 00200 {x_p_cov(l1,l2) = xpcov(l1,3+l2);} 00201 } 00202 00203 00204 // AlgebraicMatrix p_x_cov = xpcov.sub(4,7,1,3); 00205 // AlgebraicMatrix x_p_cov = xpcov.sub(1,3,4,7); 00206 00207 //here the clhep must be worken around 00208 // cout<<"p_x_cov"<< p_x_cov<<endl; 00209 // cout<<"x_p_cov"<< x_p_cov<<endl; 00210 00211 //putting everything to the joint covariance matrix: 00212 //diagonal momentum-momentum elements: 00213 cov.sub(i_int*4 + 4, i_int*4 + 4,m_m_cov); 00214 00215 //position momentum elements: 00216 cov.sub(1,i_int*4 + 4,x_p_cov); 00217 cov.sub(i_int*4 + 4,1,p_x_cov); 00218 00219 //off diagonal elements: corellations 00220 // track momentum - track momentum 00221 } 00222 int j_int = 0; 00223 for(vector<RefCountedVertexTrack>::const_iterator rt_j = refTracks.begin(); rt_j != refTracks.end(); rt_j++) 00224 { 00225 if(i_int < j_int) 00226 { 00227 AlgebraicMatrix i_k_cov_m = asHepMatrix<4,4>(vtx.tkToTkCovariance((*rt_i),(*rt_j))); 00228 // cout<<"i_k_cov_m"<<i_k_cov_m <<endl; 00229 cov.sub(i_int*4 + 4, j_int*4 + 4,i_k_cov_m); 00230 cov.sub(j_int*4 + 4, i_int*4 + 4,i_k_cov_m); 00231 } 00232 j_int++; 00233 } 00234 i_int++; 00235 } 00236 // cout<<"jac"<<jac<<endl; 00237 // cout<<"cov"<<cov<<endl; 00238 00239 return jac*cov*jac.T(); 00240 }
Definition at line 37 of file FinalTreeBuilder.h.
Referenced by buildTree(), FinalTreeBuilder(), and ~FinalTreeBuilder().
Definition at line 38 of file FinalTreeBuilder.h.
Referenced by buildTree(), FinalTreeBuilder(), and ~FinalTreeBuilder().