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