BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
VertexFit.cxx
Go to the documentation of this file.
1#include <cassert>
2
3#include "TStopwatch.h"
4#include "VertexFit/BField.h"
5#include "VertexFit/HTrackParameter.h"
6#include "VertexFit/VertexConstraints.h"
7#include "VertexFit/VertexFit.h"
8
9const int VertexFit::NTRKPAR = 6;
10const int VertexFit::NVTXPAR = 3;
11const int VertexFit::NCONSTR = 2;
12
13VertexFit* VertexFit::m_pointer = 0;
14
15VertexFit* VertexFit::instance() {
16 if ( m_pointer ) return m_pointer;
17 m_pointer = new VertexFit();
18 return m_pointer;
19}
20
21VertexFit::~VertexFit() {
22 // if (m_pointer) delete m_pointer;
23}
24
25VertexFit::VertexFit() { ; }
26
27void VertexFit::init() {
28 // derived from TrackPool
29 clearWTrackOrigin();
30 clearWTrackInfit();
31 clearWTrackList();
32 clearGammaShape();
33 clearGammaShapeList();
34 clearMapkinematic();
35 clearMappositionA();
36 clearMappositionB();
37 clearone();
38 cleartwo();
39
40 m_vpar_origin.clear();
41 m_vpar_infit.clear();
42 m_vc.clear();
43 m_chisq.clear();
44 m_chi = 9999.;
45 m_virtual_wtrk.clear();
46 m_niter = 10;
47 m_chiter = 1.0e-3;
48 m_chicut = 1000;
49
50 m_TRA = HepMatrix( 6, 7, 0 );
51 m_TRA[0][0] = 1.0;
52 m_TRA[1][1] = 1.0;
53 m_TRA[2][2] = 1.0;
54 m_TRA[3][4] = 1.0;
55 m_TRA[4][5] = 1.0;
56 m_TRA[5][6] = 1.0;
57 m_TRB = HepMatrix( 7, 6, 0 );
58 m_TRB[0][0] = 1.0;
59 m_TRB[1][1] = 1.0;
60 m_TRB[2][2] = 1.0;
61 m_TRB[4][3] = 1.0;
62 m_TRB[5][4] = 1.0;
63 m_TRB[6][5] = 1.0;
64
65 m_factor = 1.000;
66}
67
68//
69// Add Beam Fit
70//
71void VertexFit::AddBeamFit( int number, VertexParameter vpar, int n ) {
72 std::vector<int> tlis = AddList( n );
73 VertexConstraints vc;
74 vc.FixedVertexConstraints( tlis );
75 m_vc.push_back( vc );
76 m_vpar_origin.push_back( vpar );
77 m_vpar_infit.push_back( vpar );
78 if ( (unsigned int)number != m_vc.size() - 1 )
79 std::cout << "wrong kinematic constraints index" << std::endl;
80}
81
82//
83// Add Vertex
84//
85void VertexFit::AddVertex( int number, VertexParameter vpar, std::vector<int> tlis ) {
86 VertexConstraints vc;
87 vc.CommonVertexConstraints( tlis );
88 m_vc.push_back( vc );
89 HepVector vx( 3, 0 );
90 for ( unsigned int i = 0; i < tlis.size(); i++ ) vx += wTrackOrigin( tlis[i] ).X();
91 vx = vx / tlis.size();
92 VertexParameter n_vpar = vpar;
93 n_vpar.setVx( vx );
94 m_vpar_origin.push_back( n_vpar );
95 m_vpar_infit.push_back( n_vpar );
96 WTrackParameter vwtrk;
97 m_virtual_wtrk.push_back( vwtrk );
98 if ( (unsigned int)number != m_vc.size() - 1 )
99 std::cout << "wrong kinematic constraints index" << std::endl;
100}
101
102void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2 ) {
103 std::vector<int> tlis = AddList( n1, n2 );
104 AddVertex( number, vpar, tlis );
105}
106
107void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3 ) {
108 std::vector<int> tlis = AddList( n1, n2, n3 );
109 AddVertex( number, vpar, tlis );
110}
111
112void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4 ) {
113 std::vector<int> tlis = AddList( n1, n2, n3, n4 );
114 AddVertex( number, vpar, tlis );
115}
116
117void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
118 int n5 ) {
119 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5 );
120 AddVertex( number, vpar, tlis );
121}
122
123void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
124 int n5, int n6 ) {
125 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5, n6 );
126 AddVertex( number, vpar, tlis );
127}
128
129void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
130 int n5, int n6, int n7 ) {
131 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5, n6, n7 );
132 AddVertex( number, vpar, tlis );
133}
134
135void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
136 int n5, int n6, int n7, int n8 ) {
137 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5, n6, n7, n8 );
138 AddVertex( number, vpar, tlis );
139}
140
141void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
142 int n5, int n6, int n7, int n8, int n9 ) {
143 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5, n6, n7, n8, n9 );
144 AddVertex( number, vpar, tlis );
145}
146
147void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
148 int n5, int n6, int n7, int n8, int n9, int n10 ) {
149 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5, n6, n7, n8, n9, n10 );
150 AddVertex( number, vpar, tlis );
151}
152
153void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
154 int n5, int n6, int n7, int n8, int n9, int n10, int n11 ) {
155 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11 );
156 AddVertex( number, vpar, tlis );
157}
158
159void VertexFit::AddVertex( int number, VertexParameter vpar, int n1, int n2, int n3, int n4,
160 int n5, int n6, int n7, int n8, int n9, int n10, int n11,
161 int n12 ) {
162 std::vector<int> tlis = AddList( n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12 );
163 AddVertex( number, vpar, tlis );
164}
165
166bool VertexFit::Fit( int n ) {
167 bool okfit = false;
168 TStopwatch timer;
169 m_cpu = HepVector( 10, 0 );
170 if ( n < 0 || (unsigned int)n >= m_vc.size() ) return okfit;
171
172 timer.Start();
173 int ifail;
174 m_nvtrk = numberWTrack();
175
176 m_pOrigin = HepVector( m_nvtrk * NTRKPAR, 0 );
177 m_pInfit = HepVector( m_nvtrk * NTRKPAR, 0 );
178 m_pcovOrigin = HepSymMatrix( m_nvtrk * NTRKPAR, 0 );
179 m_pcovInfit = HepSymMatrix( m_nvtrk * NTRKPAR, 0 );
180
181 int ntrk = numberWTrack();
182 for ( unsigned int itk = 0; itk < ntrk; itk++ )
183 {
184 setWTrackInfit( itk, wTrackOrigin( itk ) );
185 setPOrigin( itk, Convert76( wTrackOrigin( itk ).w() ) );
186 setPCovOrigin( itk, wTrackOrigin( itk ).Ew().similarity( m_TRA ) );
187 }
188 m_pInfit = m_pOrigin;
189
190 m_xOrigin = HepVector( NVTXPAR, 0 );
191 m_xInfit = HepVector( NVTXPAR, 0 );
192 m_xcovOrigin = HepSymMatrix( NVTXPAR, 0 );
193 m_xcovInfit = HepSymMatrix( NVTXPAR, 0 );
194 m_xcovOriginInversed = HepSymMatrix( NVTXPAR, 0 );
195 m_xcovInfitInversed = HepSymMatrix( NVTXPAR, 0 );
196
197 m_xOrigin = m_vpar_origin[n].Vx();
198 m_xcovOrigin = m_vpar_origin[n].Evx();
199 m_xcovOriginInversed = m_xcovOrigin.inverse( ifail );
200 m_xInfit = m_xOrigin;
201
202 m_vpar_infit[n] = m_vpar_origin[n];
203
204 m_B = HepMatrix( NCONSTR * m_nvtrk, NVTXPAR, 0 );
205 m_A = HepMatrix( NCONSTR * m_nvtrk, NTRKPAR * m_nvtrk, 0 );
206 m_BT = HepMatrix( NVTXPAR, NCONSTR * m_nvtrk, 0 );
207 m_AT = HepMatrix( NTRKPAR * m_nvtrk, NCONSTR * m_nvtrk, 0 );
208 m_G = HepVector( NCONSTR * m_nvtrk, 0 );
209 m_W = HepSymMatrix( NCONSTR * m_nvtrk, 0 );
210 m_E = HepMatrix( NTRKPAR * m_nvtrk, NVTXPAR, 0 );
211
212 timer.Stop();
213 m_cpu[0] += timer.CpuTime();
214
215 // iteration loop
216 std::vector<double> chisq;
217 chisq.clear();
218 for ( int it = 0; it < m_niter; it++ )
219 {
220 timer.Start();
221 UpdateConstraints( m_vc[n] );
222 timer.Stop();
223 m_cpu[1] += timer.CpuTime();
224 timer.Start();
225 fitVertex( n );
226 timer.Stop();
227 m_cpu[2] += timer.CpuTime();
228 chisq.push_back( m_chisq[n] );
229 if ( it > 0 )
230 {
231 double delchi = chisq[it] - chisq[it - 1];
232 if ( fabs( delchi ) < m_chiter ) break;
233 }
234 }
235
236 /*REVISED
237 if(m_chisq[n] >= m_chicut || m_chisq[n] < 0) return okfit;
238 REVISED*/
239 if ( m_chisq[n] >= m_chicut ) return okfit;
240
241 // update vertex and its covariance
242 m_vpar_infit[n].setVx( m_xInfit );
243 m_vpar_infit[n].setEvx( m_xcovInfit );
244
245 okfit = true;
246 return okfit;
247}
248
249bool VertexFit::BeamFit( int n ) {
250 bool okfit = false;
251 if ( n < 0 || (unsigned int)n >= m_vc.size() ) return okfit;
252 for ( unsigned int i = 0; i < ( m_vc[n].Ltrk() ).size(); i++ )
253 {
254 int itk = ( m_vc[n].Ltrk() )[i];
255 setWTrackInfit( itk, wTrackOrigin( itk ) );
256 }
257 m_vpar_infit[n] = m_vpar_origin[n];
258
259 // iteration loop
260 std::vector<double> chisq;
261 chisq.clear();
262 for ( int it = 0; it < m_niter; it++ )
263 {
264 std::vector<WTrackParameter> wlis;
265 wlis.clear();
266 for ( unsigned int i = 0; i < ( m_vc[n].Ltrk() ).size(); i++ )
267 {
268 int itk = ( m_vc[n].Ltrk() )[i];
269 wlis.push_back( wTrackInfit( itk ) );
270 }
271 VertexParameter vpar = m_vpar_infit[n];
272 m_vc[n].UpdateConstraints( vpar, wlis );
273
274 fitBeam( n );
275 chisq.push_back( m_chisq[n] );
276 if ( it > 0 )
277 {
278 double delchi = chisq[it] - chisq[it - 1];
279 if ( fabs( delchi ) < m_chiter ) break;
280 }
281 }
282 if ( m_chisq[n] >= m_chicut ) return okfit;
283 swimBeam( n );
284 okfit = true;
285 return okfit;
286}
287
288bool VertexFit::Fit() {
289 bool okfit = false;
290 double mychi = 0;
291 for ( unsigned int n = 0; n < (int)( m_vc.size() ); n++ )
292 {
293 Fit( n );
294 if ( m_chisq[n] >= m_chicut ) return okfit;
295 swimVertex( n );
296 mychi = mychi + m_chisq[n];
297 }
298 m_chi = mychi;
299 okfit = true;
300 return okfit;
301}
302
303void VertexFit::fitVertex( int n ) {
304 if ( m_chisq.size() == 0 )
305 {
306 for ( unsigned int i = 0; i < m_vc.size(); i++ ) m_chisq.push_back( 9999.0 );
307 }
308 TStopwatch timer;
309 VertexConstraints vc = m_vc[n];
310
311 int ifail;
312 int NSIZE = ( vc.Ltrk() ).size();
313
314 // get new Vx
315 timer.Start();
316 m_xcovInfitInversed = m_xcovOriginInversed;
317
318 for ( unsigned int i = 0; i < NSIZE; i++ )
319 {
320 int itk = ( vc.Ltrk() )[i];
321 m_xcovInfitInversed += vfW( itk ).similarity( vfBT( itk ) );
322 }
323 m_xcovInfit = m_xcovInfitInversed.inverse( ifail );
324
325 // calculate Kq and E
326 m_KQ = HepMatrix( NVTXPAR, m_nvtrk * NCONSTR, 0 );
327 m_E = HepMatrix( m_nvtrk * NTRKPAR, NVTXPAR, 0 );
328 for ( unsigned int i = 0; i < NSIZE; i++ )
329 {
330 int itk = ( vc.Ltrk() )[i];
331 setKQ( itk, ( m_xcovInfit * vfBT( itk ) * vfW( itk ) ) );
332 setE( itk, ( -pcovOrigin( itk ) * vfAT( itk ) * vfKQ( itk ).T() ) );
333 }
334 // update vertex position
335 m_xInfit = m_xOrigin;
336 for ( unsigned int i = 0; i < NSIZE; i++ )
337 {
338 int itk = ( vc.Ltrk() )[i];
339 m_xInfit -= vfKQ( itk ) * vfG( itk );
340 }
341 // update Track parameter
342 HepVector dq0q( NVTXPAR, 0 );
343 dq0q = m_xcovInfitInversed * ( m_xOrigin - m_xInfit );
344 for ( unsigned int i = 0; i < NSIZE; i++ )
345 {
346 int itk = ( vc.Ltrk() )[i];
347 HepVector alpha( NTRKPAR, 0 );
348 alpha = pOrigin( itk ) - pcovOrigin( itk ) * vfAT( itk ) *
349 ( vfW( itk ) * vfG( itk ) - vfKQ( itk ).T() * dq0q );
350 setPInfit( itk, alpha );
351 }
352 // get chisquare value
353 m_chisq[n] = ( m_W.similarity( m_G.T() ) -
354 m_xcovInfitInversed.similarity( ( m_xInfit - m_xOrigin ).T() ) )[0][0];
355}
356
357void VertexFit::vertexCovMatrix( int n ) {
358 TStopwatch timer;
359 timer.Start();
360
361 VertexConstraints vc = m_vc[n];
362
363 unsigned int NSIZE = vc.Ltrk().size();
364 int ifail;
365 m_pcovInfit = HepSymMatrix( NTRKPAR * m_nvtrk, 0 );
366 for ( unsigned int i = 0; i < NSIZE; i++ )
367 {
368 int itk = vc.Ltrk()[i];
369 setPCovInfit( itk, pcovOrigin( itk ) -
370 vfW( itk ).similarity( pcovOrigin( itk ) * vfAT( itk ) ) );
371 }
372 m_pcovInfit = m_pcovInfit + m_xcovInfitInversed.similarity( m_E );
373
374 timer.Stop();
375 m_cpu[3] += timer.CpuTime();
376}
377
378void VertexFit::swimVertex( int n ) {
379 TStopwatch timer;
380 timer.Start();
381 // Track parameter can be expressed as:
382 //
383 // px = px0 + a*(y0 - yv)
384 // py = py0 - a*(x0 - xv)
385 // pz = pz0
386 // e = e
387 // x = xv
388 // y = yv
389 // z = zv
390 //
391 // thus p = A * a + B * vx
392 // x = vx
393 VertexConstraints vc = m_vc[n];
394 unsigned int NSIZE = vc.Ltrk().size();
395
396 HepMatrix A( 6, 6, 0 );
397 A[0][0] = 1.0;
398 A[1][1] = 1.0;
399 A[2][2] = 1.0;
400 HepMatrix B( 6, 3, 0 );
401 B[3][0] = 1.0;
402 B[4][1] = 1.0;
403 B[5][2] = 1.0;
404 HepVector w1( 6, 0 );
405 HepVector w2( 7, 0 );
406 HepSymMatrix Ew( 7, 0 );
407 HepMatrix ew1( 6, 6, 0 );
408 HepMatrix ew2( 7, 7, 0 );
409
410 for ( unsigned int i = 0; i < NSIZE; i++ )
411 {
412 int itk = vc.Ltrk()[i];
413 // double afield = VertexFitBField::instance()->getCBz(m_xInfit,
414 //pInfit(itk).sub(4, 6));
415 double afield =
416 m_factor * VertexFitBField::instance()->getCBz( m_xInfit, pInfit( itk ).sub( 4, 6 ) );
417 double a = afield * wTrackInfit( itk ).charge();
418
419 A[0][4] = a;
420 A[1][3] = -a;
421 B[0][1] = -a;
422 B[1][0] = a;
423 w1 = A * pInfit( itk ) + B * m_xInfit;
424 ew1 = pcovInfit( itk ).similarity( A ) + m_xcovInfit.similarity( B ) +
425 A * vfE( itk ) * B.T() + B * ( vfE( itk ).T() ) * A.T();
426
427 WTrackParameter wtrk = wTrackOrigin( itk );
428 w2 = Convert67( wtrk.mass(), w1 );
429 wtrk.setW( w2 );
430
431 m_TRB[3][0] = w2[0] / w2[3];
432 m_TRB[3][1] = w2[1] / w2[3];
433 m_TRB[3][2] = w2[2] / w2[3];
434
435 ew2 = m_TRB * ew1 * m_TRB.T();
436 Ew.assign( ew2 );
437 wtrk.setEw( Ew );
438 // renew parameters of input track
439 setWTrackInfit( itk, wtrk );
440 }
441 timer.Stop();
442 m_cpu[4] += timer.CpuTime();
443}
444
445bool VertexFit::pull( int n, int itk, HepVector& p ) {
446 assert( p.num_row() == 5 );
447 vertexCovMatrix( n );
448
449 WTrackParameter wtrk0, wtrk1;
450 HepVector w1( 6, 0 );
451 HepVector w2( 7, 0 );
452 HepSymMatrix ew1( 6, 0 );
453 HepSymMatrix ew2( 7, 0 );
454 wtrk0 = wTrackOrigin( itk );
455 w1 = pInfit( itk );
456 ew1 = pcovInfit( itk );
457 w2 = Convert67( wtrk0.mass(), w1 );
458 m_TRB[3][0] = w2[0] / w2[3];
459 m_TRB[3][1] = w2[1] / w2[3];
460 m_TRB[3][2] = w2[2] / w2[3];
461 ew2 = ew1.similarity( m_TRB );
462 wtrk1.setW( w2 );
463 wtrk1.setEw( ew2 );
464 wtrk1.setCharge( wtrk0.charge() );
465
466 HTrackParameter htrk0( wtrk0 );
467 HTrackParameter htrk1( wtrk1 );
468 for ( int i = 0; i < 5; i++ )
469 {
470 double del = htrk0.eHel()[i][i] - htrk1.eHel()[i][i];
471 if ( del == 0.0 ) { return false; }
472 p[i] = ( htrk0.helix()[i] - htrk1.helix()[i] ) / sqrt( abs( del ) );
473 }
474 return true;
475}
476
477void VertexFit::fitBeam( int n ) {
478 /*
479 if(m_chisq.size() == 0) {
480 for(unsigned int i = 0; i < m_vc.size(); i++)
481 m_chisq.push_back(0.0);
482 }
483
484 VertexConstraints vc = m_vc[n];
485 VertexParameter vpar = m_vpar_origin[n];
486
487 unsigned int NSIZE = vc.Ltrk().size();
488 int ifail;
489
490 // get VD, lambda
491 for(unsigned int i = 0; i < NSIZE; i++) {
492 int itk = vc.Ltrk()[i];
493 HepVector dela0(7, 0);
494 dela0 = wTrackOrigin(itk).w()-wTrackInfit(itk).w();
495 HepSymMatrix vd(2, 0);
496 vd = ((wTrackOrigin(itk).Ew()).similarity(vc.Dc()[i])).inverse(ifail);
497 HepVector lambda(2, 0);
498 lambda = vd*(vc.Dc()[i]*dela0+vc.dc()[i]);
499 vc.setVD(i, vd);
500 vc.setLambda(i, lambda);
501 }
502 // get new dela = dela0 - Va0 Dt lambda
503 // get new Va
504 // get chisq
505 HepMatrix covax(7, 3, 0);
506 m_chisq[n] = 0;
507 for(unsigned int i = 0; i < NSIZE; i++) {
508 HepVector DtL(7, 0);
509 DtL = (vc.Dc()[i]).T() * vc.lambda()[i];
510 int itk = vc.Ltrk()[i];
511 HepVector dela(7, 0);
512 HepVector dela0(7, 0);
513 dela0 = wTrackOrigin(itk).w()-wTrackInfit(itk).w();
514 WTrackParameter wtrk = wTrackOrigin(itk);
515 dela = -wTrackOrigin(itk).Ew() * DtL;
516 wtrk.setW(wtrk.w()+dela);
517 HepSymMatrix Va(7, 0);
518 HepSymMatrix Va0(7, 0);
519 Va0 = wTrackOrigin(itk).Ew();
520 HepMatrix DcT(7, 2, 0);
521 DcT = (vc.Dc()[i]).T();
522 HepSymMatrix DVdD(7, 0);
523 DVdD = (vc.VD()[i]).similarity(DcT);
524 Va = Va0 - DVdD.similarity(Va0);
525 wtrk.setEw(Va);
526 setWTrackInfit(itk, wtrk);
527 HepVector dc(2, 0);
528 dc = vc.Dc()[i]*dela0 +vc.dc()[i];
529 m_chisq[n] = m_chisq[n] + dot(vc.lambda()[i], dc);
530 }
531 m_vpar_infit[n] = vpar;
532 m_vc[n] = vc;
533 */
534 // No needed
535}
536
537void VertexFit::swimBeam( int n ) {
538 /*
539 // const double alpha = -0.00299792458;
540 //
541 // Track parameter can be expressed as:
542 //
543 // px = px0 + a*(y0 - yv)
544 // py = py0 - a*(x0 - xv)
545 // pz = pz0
546 // e = e
547 // x = xv
548 // y = yv
549 // z = zv
550 //
551 // thus p = A * a + B * vx
552 // x = vx
553
554
555 VertexConstraints vc = m_vc[n];
556 VertexParameter vpar = m_vpar_infit[n];
557 unsigned int NSIZE = vc.Ltrk().size();
558
559 for(unsigned int i = 0; i < NSIZE; i++) {
560 int itk = vc.Ltrk()[i];
561 HepMatrix A(4, 7, 0);
562 HepMatrix B(4, 3, 0);
563
564 double afield = VertexFitBField::instance()->getCBz(vpar.Vx(), pInfit(itk).sub(4, 6));
565 double a = afield * wTrackInfit(itk).charge();
566 A[0][0] = 1.0;
567 A[0][5] = a;
568 A[1][1] = 1.0;
569 A[1][4] = -a;
570 A[2][2] = 1.0;
571 A[3][3] = 1.0;
572 B[0][1] = -a;
573 B[1][0] = a;
574 HepVector p(4, 0);
575 p = A * wTrackInfit(itk).w() + B * vpar.Vx();
576 HepVector x(3, 0);
577 x = vpar.Vx();
578 HepVector w(7, 0);
579 HepSymMatrix Ew(7, 0);
580 for(int j = 0; j < 4; j++)
581 w[j] = p[j];
582 for(int j = 0; j < 3; j++)
583 w[j+4] = x[j];
584
585 WTrackParameter wtrk;
586 wtrk.setCharge(wTrackInfit(itk).charge());
587 wtrk.setW(w);
588 HepSymMatrix Vpv(4, 0);
589 Vpv = (wTrackInfit(itk).Ew()).similarity(A);
590 for(int j = 0; j < 4; j++)
591 for(int k= 0; k < 4; k++)
592 Ew[j][k] = Vpv[j][k];
593 wtrk.setEw(Ew);
594 setWTrackInfit(itk, wtrk);
595 }
596 */
597 // No needed
598}
599
600void VertexFit::BuildVirtualParticle( int n ) {
601
602 vertexCovMatrix( n );
603 TStopwatch timer;
604 timer.Start();
605
606 HepMatrix A( NTRKPAR, NTRKPAR * m_nvtrk, 0 );
607 HepMatrix B( NTRKPAR, NVTXPAR, 0 );
608 VertexConstraints vc = m_vc[n];
609 unsigned int NSIZE = vc.Ltrk().size();
610 int charge = 0;
611
612 HepMatrix Ai( 6, 6, 0 );
613 Ai[0][0] = 1.0;
614 Ai[1][1] = 1.0;
615 Ai[2][2] = 1.0;
616 HepMatrix Bi( 6, 3, 0 );
617 Bi[3][0] = 1.0;
618 Bi[4][1] = 1.0;
619 Bi[5][2] = 1.0;
620 HepVector w1( 6, 0 );
621 HepVector w2( 7, 0 );
622 HepSymMatrix Ew( 7, 0 );
623 HepMatrix ew1( 6, 6, 0 );
624 HepMatrix ew2( 7, 7, 0 );
625 double totalE = 0;
626
627 for ( unsigned int i = 0; i < NSIZE; i++ )
628 {
629 int itk = vc.Ltrk()[i];
630 charge += wTrackInfit( itk ).charge();
631 // double afield = VertexFitBField::instance()->getCBz(m_xInfit,
632 //pInfit(itk).sub(4, 6));
633 double afield =
634 m_factor * VertexFitBField::instance()->getCBz( m_xInfit, pInfit( itk ).sub( 4, 6 ) );
635 double a = afield * wTrackInfit( itk ).charge();
636
637 // totalE += wTrackOrigin(itk).w()[3];
638 totalE += wTrackInfit( itk ).w()[3]; // update total energy after fit, commit by sunhk,
639 // thanks to Chengzhi He.
640 Ai[0][4] = a;
641 Ai[1][3] = -a;
642 Bi[0][1] = -a;
643 Bi[1][0] = a;
644 A.sub( 1, NTRKPAR * itk + 1, Ai );
645 B += Bi;
646 }
647 B[3][0] = 1.0;
648 B[4][1] = 1.0;
649 B[5][2] = 1.0;
650
651 w1 = A * m_pInfit + B * m_xInfit;
652 ew1 = m_pcovInfit.similarity( A ) + m_xcovInfit.similarity( B ) + A * m_E * B.T() +
653 B * ( m_E.T() ) * A.T();
654
655 // convert w1(6x1) to w2(7x1)
656 w2[0] = w1[0];
657 w2[1] = w1[1];
658 w2[2] = w1[2];
659 w2[3] = totalE;
660 w2[4] = w1[3];
661 w2[5] = w1[4];
662 w2[6] = w1[5];
663 // convert ew1(6x6) to ew2(7x7)
664 m_TRB[3][0] = w1[0] / totalE;
665 m_TRB[3][1] = w1[1] / totalE;
666 m_TRB[3][2] = w1[2] / totalE;
667 ew2 = m_TRB * ew1 * m_TRB.T();
668 Ew.assign( ew2 );
669 WTrackParameter vwtrk;
670 vwtrk.setCharge( charge );
671 vwtrk.setW( w2 );
672 vwtrk.setEw( Ew );
673
674 m_virtual_wtrk[n] = vwtrk;
675 timer.Stop();
676 m_cpu[5] += timer.CpuTime();
677}
678
679void VertexFit::UpdateConstraints( const VertexConstraints& vc ) {
680 int ntrk = ( vc.Ltrk() ).size();
681 int type = vc.type();
682 switch ( type )
683 {
684 case 1: {
685 for ( unsigned int i = 0; i < ntrk; i++ )
686 {
687 int itk = ( vc.Ltrk() )[i];
688 HepVector alpha( 6, 0 );
689 double mass, e;
690 HepLorentzVector p;
691 HepPoint3D x, vx;
692 alpha = pInfit( itk );
693
694 mass = wTrackOrigin( itk ).mass();
695 e = sqrt( mass * mass + alpha[0] * alpha[0] + alpha[1] * alpha[1] +
696 alpha[2] * alpha[2] );
697 p = HepLorentzVector( alpha[0], alpha[1], alpha[2], e );
698 x = HepPoint3D( alpha[3], alpha[4], alpha[5] );
699
700 vx = HepPoint3D( m_xInfit[0], m_xInfit[1], m_xInfit[2] );
701 HepPoint3D delx = vx - x;
702
703 // double afield =
704 //VertexFitBField::instance()->getCBz(m_xInfit, wTrackOrigin(itk).X());
705 double afield =
706 m_factor * VertexFitBField::instance()->getCBz( m_xInfit, wTrackOrigin( itk ).X() );
707 double a = afield * ( 0.0 + wTrackOrigin( itk ).charge() );
708
709 double J = a * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
710 J = std::min( J, 1 - 1e-4 );
711 J = std::max( J, -1 + 1e-4 );
712 double Rx =
713 delx.x() - 2 * p.px() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
714 double Ry =
715 delx.y() - 2 * p.py() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
716 double S = 1.0 / sqrt( 1 - J * J ) / p.perp2();
717 // dc
718 HepVector dc( 2, 0 );
719 dc[0] = delx.y() * p.px() - delx.x() * p.py() -
720 0.5 * a * ( delx.x() * delx.x() + delx.y() * delx.y() );
721 dc[1] = delx.z() - p.pz() / a * asin( J );
722 // setd(itk, dc);
723 // Ec
724 HepMatrix Ec( 2, 3, 0 );
725 // m_Ec.push_back(Ec);
726 // Dc
727 HepMatrix Dc( 2, 6, 0 );
728 Dc[0][0] = delx.y();
729 Dc[0][1] = -delx.x();
730 Dc[0][2] = 0;
731 Dc[0][3] = p.py() + a * delx.x();
732 Dc[0][4] = -p.px() + a * delx.y();
733 Dc[0][5] = 0;
734 Dc[1][0] = -p.pz() * S * Rx;
735 Dc[1][1] = -p.pz() * S * Ry;
736 Dc[1][2] = -asin( J ) / a;
737 Dc[1][3] = p.px() * p.pz() * S;
738 Dc[1][4] = p.py() * p.pz() * S;
739 Dc[1][5] = -1.0;
740 // setD(itk, Dc);
741 // setDT(itk, Dc.T());
742 // VD
743 HepSymMatrix vd( 2, 0 );
744 int ifail;
745 vd = pcovOrigin( itk ).similarity( Dc );
746 vd.invert( ifail );
747 // setVD(itk, vd);
748 }
749 break;
750 }
751 case 2:
752 default: {
753 for ( unsigned int i = 0; i < ntrk; i++ )
754 {
755 int itk = ( vc.Ltrk() )[i];
756 HepVector alpha( 6, 0 );
757 double mass, e;
758 HepLorentzVector p;
759 HepPoint3D x, vx;
760 alpha = pInfit( itk );
761 // p = HepLorentzVector(alpha[0], alpha[1], alpha[2], alpha[3]);
762 // x = HepPoint3D(alpha[4], alpha[5], alpha[6]);
763
764 mass = wTrackOrigin( itk ).mass();
765 e = sqrt( mass * mass + alpha[0] * alpha[0] + alpha[1] * alpha[1] +
766 alpha[2] * alpha[2] );
767 p = HepLorentzVector( alpha[0], alpha[1], alpha[2], e );
768 x = HepPoint3D( alpha[3], alpha[4], alpha[5] );
769
770 vx = HepPoint3D( m_xInfit[0], m_xInfit[1], m_xInfit[2] );
771 HepPoint3D delx = vx - x;
772
773 if ( wTrackOrigin( itk ).charge() == 0 )
774 {
775 // dc -->g
776 HepVector dc( 2, 0 );
777 dc[0] = p.pz() * delx.x() - p.px() * delx.z();
778 dc[1] = p.pz() * delx.y() - p.py() * delx.z();
779 // Ec --> B
780 HepMatrix Ec( 2, 3, 0 );
781 Ec[0][0] = p.pz();
782 Ec[0][1] = 0;
783 Ec[0][2] = -p.px();
784 Ec[1][0] = 0;
785 Ec[1][1] = p.pz();
786 Ec[1][2] = -p.py();
787 setB( itk, Ec );
788 setBT( itk, Ec.T() );
789 // Dc
790 HepMatrix Dc( 2, 6, 0 );
791 Dc[0][0] = -delx.z();
792 Dc[0][1] = 0;
793 Dc[0][2] = delx.x();
794 Dc[0][3] = -p.pz();
795 Dc[0][4] = 0;
796 Dc[0][5] = p.px();
797 Dc[1][0] = 0;
798 Dc[1][1] = -delx.z();
799 Dc[1][2] = delx.y();
800 Dc[1][3] = 0;
801 Dc[1][4] = -p.pz();
802 Dc[1][5] = p.py();
803 setA( itk, Dc );
804 setAT( itk, Dc.T() );
805
806 // VD --> W
807 HepSymMatrix vd( 2, 0 );
808 int ifail;
809
810 vd = pcovOrigin( itk ).similarity( Dc );
811 vd.invert( ifail );
812 setW( itk, vd );
813
814 // G=A(p0-pe)+B(x0-xe)+d
815 HepVector gc( 2, 0 );
816 gc = Dc * ( pOrigin( itk ) - pInfit( itk ) ) + Ec * ( m_xOrigin - m_xInfit ) + dc;
817 setG( itk, gc );
818 }
819 else
820 {
821 // double afield =
822 //VertexFitBField::instance()->getCBz(m_xInfit,wTrackOrigin(itk).X());
823 double afield = m_factor * VertexFitBField::instance()->getCBz(
824 m_xInfit, wTrackOrigin( itk ).X() );
825 double a = afield * ( 0.0 + wTrackOrigin( itk ).charge() );
826 double J = a * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
827 J = std::min( J, 1 - 1e-4 );
828 J = std::max( J, -1 + 1e-4 );
829 double Rx =
830 delx.x() - 2 * p.px() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
831 double Ry =
832 delx.y() - 2 * p.py() * ( delx.x() * p.px() + delx.y() * p.py() ) / p.perp2();
833 double S = 1.0 / sqrt( 1 - J * J ) / p.perp2();
834
835 // dc
836 HepVector dc( 2, 0 );
837 dc[0] = delx.y() * p.px() - delx.x() * p.py() -
838 0.5 * a * ( delx.x() * delx.x() + delx.y() * delx.y() );
839 dc[1] = delx.z() - p.pz() / a * asin( J );
840 // Ec
841 HepMatrix Ec( 2, 3, 0 );
842 Ec[0][0] = -p.py() - a * delx.x();
843 Ec[0][1] = p.px() - a * delx.y();
844 Ec[0][2] = 0;
845 Ec[1][0] = -p.px() * p.pz() * S;
846 Ec[1][1] = -p.py() * p.pz() * S;
847 Ec[1][2] = 1.0;
848 setB( itk, Ec );
849 setBT( itk, Ec.T() );
850
851 // Dc
852 HepMatrix Dc( 2, 6, 0 );
853 Dc[0][0] = delx.y();
854 Dc[0][1] = -delx.x();
855 Dc[0][2] = 0;
856 Dc[0][3] = p.py() + a * delx.x();
857 Dc[0][4] = -p.px() + a * delx.y();
858 Dc[0][5] = 0;
859
860 Dc[1][0] = -p.pz() * S * Rx;
861 Dc[1][1] = -p.pz() * S * Ry;
862 Dc[1][2] = -asin( J ) / a;
863 Dc[1][3] = p.px() * p.pz() * S;
864 Dc[1][4] = p.py() * p.pz() * S;
865 Dc[1][5] = -1.0;
866 setA( itk, Dc );
867 setAT( itk, Dc.T() );
868 // VD
869 HepSymMatrix vd( 2, 0 );
870 int ifail;
871 vd = pcovOrigin( itk ).similarity( Dc );
872 vd.invert( ifail );
873 setW( itk, vd );
874 // G = A(p0-pe)+B(x0-xe)+d
875 HepVector gc( 2, 0 );
876 gc = Dc * ( pOrigin( itk ) - pInfit( itk ) ) + Ec * ( m_xOrigin - m_xInfit ) + dc;
877 setG( itk, gc );
878 }
879 }
880 break;
881 }
882 }
883}
884
885HepVector VertexFit::Convert67( const double& mass, const HepVector& p ) {
886 // assert(p.num_row() == 6);
887 HepVector m( 7, 0 );
888 m.sub( 1, p.sub( 1, 3 ) );
889 m.sub( 5, p.sub( 4, 6 ) );
890 m[3] = sqrt( mass * mass + p[0] * p[0] + p[1] * p[1] + p[2] * p[2] );
891 return m;
892}
893
894HepVector VertexFit::Convert76( const HepVector& p ) {
895 // assert(p.num_row() == 7);
896 HepVector m( 6, 0 );
897 m.sub( 1, p.sub( 1, 3 ) );
898 m.sub( 4, p.sub( 5, 7 ) );
899 return m;
900}
mass
double mass
Definition CosmicGenerator.cxx:128
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition CosmicGenerator.h:36
n
const Int_t n
Definition DataBase/tau_mode.c:57
x
Double_t x[10]
Definition DataBase/tau_mode.c:49
alpha
double alpha
Definition EvtPhokharaDef.hh:17
w
double w
Definition EvtPhokharaDef.hh:28
dc
int dc[18]
Definition EvtPycont.cc:63
n2
int n2
Definition SD0Tag.cxx:59
n1
int n1
Definition SD0Tag.cxx:58
TrackPool::AddList
std::vector< int > AddList(int n1)
Definition TrackPool.cxx:441
TrackPool::wTrackInfit
std::vector< WTrackParameter > wTrackInfit() const
Definition Analysis/VertexFit/include/VertexFit/TrackPool.h:72
TrackPool::wTrackOrigin
std::vector< WTrackParameter > wTrackOrigin() const
Definition Analysis/VertexFit/include/VertexFit/TrackPool.h:71
TrackPool::setWTrackInfit
void setWTrackInfit(const int n, const WTrackParameter wtrk)
Definition Analysis/VertexFit/include/VertexFit/TrackPool.h:105
VertexConstraints::FixedVertexConstraints
void FixedVertexConstraints(std::vector< int > tlis)
Definition VertexConstraints.cxx:27
VertexConstraints::CommonVertexConstraints
void CommonVertexConstraints(std::vector< int > tlis)
Definition VertexConstraints.cxx:22
VertexFitBField::getCBz
double getCBz(const HepVector &vtx, const HepVector &trackPosition)
Definition Analysis/VertexFit/src/BField.cxx:65
VertexFit::wtrk
WTrackParameter wtrk(int n) const
Definition Analysis/VertexFit/include/VertexFit/VertexFit.h:91
VertexFit::BeamFit
bool BeamFit(int n)
Definition VertexFit.cxx:249
VertexFit::AddBeamFit
void AddBeamFit(int number, VertexParameter vpar, int n)
Definition VertexFit.cxx:71
VertexFit::init
void init()
Definition VertexFit.cxx:27
VertexFit::AddVertex
void AddVertex(int number, VertexParameter vpar, std::vector< int > lis)
Definition VertexFit.cxx:85
VertexFit::pull
bool pull(int n, int itk, HepVector &p)
Definition VertexFit.cxx:445
VertexFit::instance
static VertexFit * instance()
Definition VertexFit.cxx:15
VertexFit::BuildVirtualParticle
void BuildVirtualParticle(int number)
Definition VertexFit.cxx:600
VertexFit::Fit
bool Fit()
Definition VertexFit.cxx:288