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