BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
VertexConstraints.cxx
Go to the documentation of this file.
1#include "VertexFit/VertexConstraints.h"
2#include "VertexFit/BField.h"
3#include "VertexFit/TrackPool.h"
4
5const double VertexConstraints ::Alpha = -0.00299792458;
6const int VertexConstraints::FixedVertex = 1;
7const int VertexConstraints::CommonVertex = 2;
8
10 m_Ec.clear();
11 m_Dc.clear();
12 m_dc.clear();
13 m_VD.clear();
14 m_EVDE.clear();
15 m_lambda0.clear();
16 m_lambda.clear();
17 m_covax.clear();
18 m_ltrk.clear();
19 m_type = 0;
20}
21
22void VertexConstraints::CommonVertexConstraints( std::vector<int> tlis ) {
23 m_ltrk = tlis;
25}
26
27void VertexConstraints::FixedVertexConstraints( std::vector<int> tlis ) {
28 m_ltrk = tlis;
30}
31
33 std::vector<WTrackParameter> wlis ) {
34 m_Ec.clear();
35 m_Dc.clear();
36 m_dc.clear();
37 m_VD.clear();
38 m_EVDE.clear();
39 m_lambda0.clear();
40 m_lambda.clear();
41 m_covax.clear();
42
43 switch ( m_type )
44 {
45 case FixedVertex: {
46 WTrackParameter wtrk;
47 for ( unsigned int i = 0; i < wlis.size(); i++ )
48 {
49 wtrk = wlis[i];
50 HepLorentzVector p = wtrk.p();
51 HepPoint3D x = wtrk.x();
52 HepPoint3D delx = vpar.vx() - x;
53
54 double afield = VertexFitBField::instance()->getCBz( vpar.Vx(), wtrk.X() );
55 double a = afield * ( 0.0 + wtrk.charge() );
56 double J = a * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
57
58 J = std::min( J, 1 - 1e-4 );
59 J = std::max( J, -1 + 1e-4 );
60 double Rx =
61 delx.x() - 2 * p.px() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
62 double Ry =
63 delx.y() - 2 * p.py() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
64 double S = 1.0 / sqrt( 1 - J * J ) / p.perp2();
65
66 // dc
67 HepVector dc( 2, 0 );
68 dc[0] = delx.y() * p.px() - delx.x() * p.py() -
69 0.5 * a * ( delx.x() * delx.x() + delx.y() * delx.y() );
70 dc[1] = delx.z() - p.pz() / a * asin( J );
71 m_dc.push_back( dc );
72 // Ec
73 HepMatrix Ec( 2, 3, 0 );
74 m_Ec.push_back( Ec );
75 // Dc
76 HepMatrix Dc( 2, 6, 0 );
77 Dc[0][0] = delx.y();
78 Dc[0][1] = -delx.x();
79 Dc[0][2] = 0;
80 Dc[0][3] = p.py() + a * delx.x();
81 Dc[0][4] = -p.px() + a * delx.y();
82 Dc[0][5] = 0;
83 Dc[1][0] = -p.pz() * S * Rx;
84 Dc[1][1] = -p.pz() * S * Ry;
85 Dc[1][2] = -asin( J ) / a;
86 Dc[1][3] = p.px() * p.pz() * S;
87 Dc[1][4] = p.py() * p.pz() * S;
88 Dc[1][5] = -1.0;
89 m_Dc.push_back( Dc );
90 // VD
91 HepSymMatrix VD( 2, 0 );
92 // int ifail;
93 // VD = ((wtrk.Ew()).similarity(Dc)).inverse(ifail); // 2x2 matrix (Dc Ew Dc.T())
94 m_VD.push_back( VD );
95 // EVDE
96 HepSymMatrix EVDE( 3, 0 );
97 m_EVDE.push_back( EVDE );
98 // lambda0
99 HepVector lambda0( 2, 0 );
100 m_lambda0.push_back( lambda0 );
101 // lambda
102 HepVector lambda( 2, 0 );
103 m_lambda.push_back( lambda );
104 // covax
105 HepMatrix covax( 6, 3, 0 );
106 m_covax.push_back( covax );
107 break;
108 }
109 }
110 case CommonVertex:
111 default: {
112 WTrackParameter wtrk;
113 for ( unsigned int i = 0; i < wlis.size(); i++ )
114 {
115 wtrk = wlis[i];
116 HepLorentzVector p = wtrk.p();
117 HepPoint3D x = wtrk.x();
118 HepPoint3D delx = vpar.vx() - x;
119 if ( wtrk.charge() == 0 )
120 {
121 // dc
122 HepVector dc( 2, 0 );
123 dc[0] = p.pz() * delx.x() - p.px() * delx.z();
124 dc[1] = p.pz() * delx.y() - p.py() * delx.z();
125 m_dc.push_back( dc );
126 // Ec
127 HepMatrix Ec( 2, 3, 0 );
128 Ec[0][0] = p.pz();
129 Ec[0][1] = 0;
130 Ec[0][2] = -p.px();
131 Ec[1][0] = 0;
132 Ec[1][1] = p.pz();
133 Ec[1][2] = -p.py();
134 m_Ec.push_back( Ec );
135 // Dc
136 HepMatrix Dc( 2, 6, 0 );
137 Dc[0][0] = -delx.z();
138 Dc[0][1] = 0;
139 Dc[0][2] = delx.x();
140 Dc[0][3] = -p.pz();
141 Dc[0][4] = 0;
142 Dc[0][5] = p.px();
143 Dc[1][0] = 0;
144 Dc[1][1] = -delx.z();
145 Dc[1][2] = delx.y();
146 Dc[1][3] = 0;
147 Dc[1][4] = -p.pz();
148 Dc[1][5] = p.py();
149 m_Dc.push_back( Dc );
150 }
151 else
152 {
153 double afield = VertexFitBField::instance()->getCBz( vpar.Vx(), wtrk.X() );
154 double a = afield * ( 0.0 + wtrk.charge() );
155 double J = a * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
156 J = std::min( J, 1 - 1e-4 );
157 J = std::max( J, -1 + 1e-4 );
158 double Rx =
159 delx.x() - 2 * p.px() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
160 double Ry =
161 delx.y() - 2 * p.py() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
162 double S = 1.0 / sqrt( 1 - J * J ) / p.perp2();
163 // dc
164 HepVector dc( 2, 0 );
165 dc[0] = delx.y() * p.px() - delx.x() * p.py() -
166 0.5 * a * ( delx.x() * delx.x() + delx.y() * delx.y() );
167 dc[1] = delx.z() - p.pz() / a * asin( J );
168 m_dc.push_back( dc );
169 // Ec
170 HepMatrix Ec( 2, 3, 0 );
171 Ec[0][0] = -p.py() - a * delx.x();
172 Ec[0][1] = p.px() - a * delx.y();
173 Ec[0][2] = 0;
174 Ec[1][0] = -p.px() * p.pz() * S;
175 Ec[1][1] = -p.py() * p.pz() * S;
176 Ec[1][2] = 1.0;
177 m_Ec.push_back( Ec );
178 // Dc
179 HepMatrix Dc( 2, 6, 0 );
180 Dc[0][0] = delx.y();
181 Dc[0][1] = -delx.x();
182 Dc[0][2] = 0;
183 Dc[0][3] = p.py() + a * delx.x();
184 Dc[0][4] = -p.px() + a * delx.y();
185 Dc[0][5] = 0;
186 Dc[1][0] = -p.pz() * S * Rx;
187 Dc[1][1] = -p.pz() * S * Ry;
188 Dc[1][2] = -asin( J ) / a;
189 Dc[1][3] = p.px() * p.pz() * S;
190 Dc[1][4] = p.py() * p.pz() * S;
191 Dc[1][5] = -1.0;
192 m_Dc.push_back( Dc );
193 }
194 // VD
195 HepSymMatrix VD( 2, 0 );
196 m_VD.push_back( VD );
197 // EVDE
198 HepSymMatrix EVDE( 3, 0 );
199 m_EVDE.push_back( EVDE );
200 // lambda0
201 HepVector lambda0( 2, 0 );
202 m_lambda0.push_back( lambda0 );
203 // lambda
204 HepVector lambda( 2, 0 );
205 m_lambda.push_back( lambda );
206 // covax
207 HepMatrix covax( 6, 3, 0 );
208 m_covax.push_back( covax );
209 }
210 break;
211 }
212 }
213}
void FixedVertexConstraints(std::vector< int > tlis)
void CommonVertexConstraints(std::vector< int > tlis)
void UpdateConstraints(VertexParameter vpar, std::vector< WTrackParameter > wlis)
double getCBz(const HepVector &vtx, const HepVector &trackPosition)