BOSS 6.6.4.p01
BESIII Offline Software System
Loading...
Searching...
No Matches
KinematicFit.cxx
Go to the documentation of this file.
1#include "VertexFit/KinematicFit.h"
2#include "VertexFit/KinematicConstraints.h"
3#include "VertexFit/HTrackParameter.h"
4#include "TStopwatch.h"
5
6
7const int KinematicFit::NTRKPAR = 3; //modify by yanl 2010.7.26 4---->3
8
9const int KinematicFit::Resonance = 1;
10const int KinematicFit::TotalEnergy = 2;
11const int KinematicFit::TotalMomentum = 4;
12const int KinematicFit::Position = 8;
13const int KinematicFit::ThreeMomentum = 16;
14const int KinematicFit::FourMomentum = 32;
15const int KinematicFit::EqualMass = 64;
16
17
18KinematicFit * KinematicFit::m_pointer = 0;
19
20KinematicFit * KinematicFit::instance() {
21 if(m_pointer) return m_pointer;
22 m_pointer = new KinematicFit();
23 return m_pointer;
24}
25
26KinematicFit::KinematicFit(){;}
27
28KinematicFit::~KinematicFit() {
29 // if(m_pointer) delete m_pointer;
30}
31
32
33void KinematicFit::init() {
34 clearWTrackOrigin();
35 clearWTrackInfit();
36 clearWTrackList();
37 //For Virtual Particles
38 //gamma shape
39 clearGammaShape();
40 clearGammaShapeList();
41 clearMapkinematic();
42 clearMappositionA();
43 clearMappositionB();
44 clearone();
45 cleartwo();
46 setBeamPosition(HepPoint3D(0.0,0.0,0.0));
47 setVBeamPosition(HepSymMatrix(3,0));
48
49 //=============
50 m_kc.clear();
51 m_chisq.clear();
52 m_chi = 9999.;
53 m_niter = 10;
54 m_chicut = 200.;
55 m_chiter = 0.005;
56 m_espread = 0.0;
57 m_kalman = 0;
58 m_collideangle = 11e-3;
59 m_flag = 0;
60 m_dynamicerror = 0;
61 m_nc = 0;
62 m_cpu = HepVector(10, 0);
63 m_massvector = HepVector(12,0);
64 m_virtual_wtrk.clear();
65}
66
67
68//
69// Add Resonance
70//
71
72void KinematicFit::AddResonance(int number, double mres, int n1) {
73 std::vector<int> tlis = AddList(n1);
74 AddResonance(number, mres, tlis);
75}
76
77void KinematicFit::AddResonance(int number, double mres, int n1, int n2) {
78 std::vector<int> tlis = AddList(n1, n2);
79 AddResonance(number, mres, tlis);
80}
81
82void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3) {
83 std::vector<int> tlis = AddList(n1, n2, n3);
84 AddResonance(number, mres, tlis);
85}
86
87void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4) {
88 std::vector<int> tlis = AddList(n1, n2, n3, n4);
89 AddResonance(number, mres, tlis);
90}
91
92void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
93 int n5) {
94 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5);
95 AddResonance(number, mres, tlis);
96}
97
98void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
99 int n5, int n6) {
100 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6);
101 AddResonance(number, mres, tlis);
102}
103
104void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
105 int n5, int n6, int n7) {
106 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7);
107 AddResonance(number, mres, tlis);
108}
109
110void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
111 int n5, int n6, int n7, int n8) {
112 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8);
113 AddResonance(number, mres, tlis);
114}
115
116void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
117 int n5, int n6, int n7, int n8, int n9) {
118 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9);
119 AddResonance(number, mres, tlis);
120}
121
122void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
123 int n5, int n6, int n7, int n8, int n9, int n10) {
124 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10);
125 AddResonance(number, mres, tlis);
126}
127
128
129void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
130 int n5, int n6, int n7, int n8, int n9,
131 int n10, int n11) {
132 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11);
133 AddResonance(number, mres, tlis);
134}
135
136void KinematicFit::AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
137 int n5, int n6, int n7, int n8, int n9,
138 int n10, int n11, int n12) {
139 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12);
140 AddResonance(number, mres, tlis);
141}
142
143void KinematicFit::AddResonance(int number, double mres, std::vector<int> tlis) {
144 KinematicConstraints kc;
145 HepSymMatrix Vre = HepSymMatrix(1,0);
146 kc.ResonanceConstraints(mres, tlis, Vre);
147 m_kc.push_back(kc);
148 if((unsigned int) number != m_kc.size()-1)
149 std::cout << "wrong kinematic constraints index" << std::endl;
150}
151
152
153//
154// Add TotalMomentum
155//
156
157void KinematicFit::AddTotalMomentum(int number, double ptot, int n1) {
158 std::vector<int> tlis = AddList(n1);
159 AddTotalMomentum(number, ptot, tlis);
160}
161void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2) {
162 std::vector<int> tlis = AddList(n1, n2);
163 AddTotalMomentum(number, ptot, tlis);
164}
165
166void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3) {
167 std::vector<int> tlis = AddList(n1, n2, n3);
168 AddTotalMomentum(number, ptot, tlis);
169}
170
171void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4) {
172 std::vector<int> tlis = AddList(n1, n2, n3, n4);
173 AddTotalMomentum(number, ptot, tlis);
174}
175
176void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
177 int n5) {
178 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5);
179 AddTotalMomentum(number, ptot, tlis);
180}
181
182void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
183 int n5, int n6) {
184 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6);
185 AddTotalMomentum(number, ptot, tlis);
186}
187
188void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
189 int n5, int n6, int n7) {
190 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7);
191 AddTotalMomentum(number, ptot, tlis);
192}
193
194void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
195 int n5, int n6, int n7, int n8) {
196 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8);
197 AddTotalMomentum(number, ptot, tlis);
198}
199
200void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
201 int n5, int n6, int n7, int n8, int n9) {
202 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9);
203 AddTotalMomentum(number, ptot, tlis);
204}
205
206void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
207 int n5, int n6, int n7, int n8, int n9,
208 int n10) {
209 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10);
210 AddTotalMomentum(number, ptot, tlis);
211}
212
213
214void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
215 int n5, int n6, int n7, int n8, int n9,
216 int n10, int n11) {
217 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11);
218 AddTotalMomentum(number, ptot, tlis);
219}
220
221void KinematicFit::AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
222 int n5, int n6, int n7, int n8, int n9,
223 int n10, int n11, int n12) {
224 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12);
225 AddTotalMomentum(number, ptot, tlis);
226}
227
228void KinematicFit::AddTotalMomentum(int number, double ptot, std::vector<int> tlis) {
229 KinematicConstraints kc;
230 kc.TotalMomentumConstraints(ptot, tlis);
231 m_kc.push_back(kc);
232 if((unsigned int) number != m_kc.size()-1)
233 std::cout << "wrong kinematic constraints index" << std::endl;
234}
235
236//
237// Add TotalEnergy
238//
239
240void KinematicFit::AddTotalEnergy(int number, double etot, int n1) {
241 std::vector<int> tlis = AddList(n1);
242 AddTotalEnergy(number, etot, tlis);
243}
244void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2) {
245 std::vector<int> tlis = AddList(n1, n2);
246 AddTotalEnergy(number, etot, tlis);
247}
248
249void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3) {
250 std::vector<int> tlis = AddList(n1, n2, n3);
251 AddTotalEnergy(number, etot, tlis);
252}
253
254void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4) {
255 std::vector<int> tlis = AddList(n1, n2, n3, n4);
256 AddTotalEnergy(number, etot, tlis);
257}
258
259void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
260 int n5) {
261 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5);
262 AddTotalEnergy(number, etot, tlis);
263}
264
265void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
266 int n5, int n6) {
267 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6);
268 AddTotalEnergy(number, etot, tlis);
269}
270
271void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
272 int n5, int n6, int n7) {
273 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7);
274 AddTotalEnergy(number, etot, tlis);
275}
276
277void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
278 int n5, int n6, int n7, int n8) {
279 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8);
280 AddTotalEnergy(number, etot, tlis);
281}
282
283void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
284 int n5, int n6, int n7, int n8, int n9) {
285 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9);
286 AddTotalEnergy(number, etot, tlis);
287}
288
289void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
290 int n5, int n6, int n7, int n8, int n9,
291 int n10) {
292 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10);
293 AddTotalEnergy(number, etot, tlis);
294}
295
296
297void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
298 int n5, int n6, int n7, int n8, int n9,
299 int n10, int n11) {
300 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11);
301 AddTotalEnergy(number, etot, tlis);
302}
303
304void KinematicFit::AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
305 int n5, int n6, int n7, int n8, int n9,
306 int n10, int n11, int n12) {
307 std::vector<int> tlis = AddList(n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11, n12);
308 AddTotalEnergy(number, etot, tlis);
309}
310
311void KinematicFit::AddTotalEnergy(int number, double etot, std::vector<int> tlis) {
312 KinematicConstraints kc;
313 kc.TotalEnergyConstraints(etot, tlis);
314 m_kc.push_back(kc);
315 if((unsigned int) number != m_kc.size()-1)
316 std::cout << "wrong kinematic constraints index" << std::endl;
317}
318//
319// Equal Mass constraints
320//
321
322void KinematicFit::AddEqualMass(int number, std::vector<int> tlis1, std::vector<int> tlis2) {
323 KinematicConstraints kc;
324 HepSymMatrix Vne = HepSymMatrix(1,0);
325 kc.EqualMassConstraints(tlis1, tlis2, Vne);
326 m_kc.push_back(kc);
327 if((unsigned int) number != m_kc.size()-1)
328 std::cout << "wrong kinematic constraints index" << std::endl;
329}
330
331//
332// Total 3-momentum constraints
333//
334
335void KinematicFit::AddThreeMomentum(int number, Hep3Vector p3) {
336 std::vector<int> tlis;
337 tlis.clear();
338 WTrackParameter wtrk;
339 KinematicConstraints kc;
340
341 for(int i = 0; i < numberWTrack(); i++) {
342 tlis.push_back(i);
343 }
344 kc.ThreeMomentumConstraints(p3, tlis);
345 m_kc.push_back(kc);
346 if((unsigned int) number != m_kc.size()-1)
347 std::cout << "wrong kinematic constraints index" << std::endl;
348}
349
350//
351// Total 4-momentum constraints
352//
353
354void KinematicFit::AddFourMomentum(int number, HepLorentzVector p4) {
355
356 std::vector<int> tlis;
357 tlis.clear();
358 KinematicConstraints kc;
359
360 for(int i = 0; i < numberWTrack(); i++) {
361 tlis.push_back(i);
362 }
363 // for(int i = 0; i < numberWTrack_V(); i++) {
364 // tlis_V.push_back(i);
365 // }
366
367 HepSymMatrix Vme = HepSymMatrix(4,0);
368 Vme[0][0] = 2*m_espread*m_espread*sin(m_collideangle)*sin(m_collideangle);
369 Vme[0][3] = 2*m_espread*m_espread*sin(m_collideangle);
370 Vme[3][3] = 2*m_espread*m_espread;
371
372 // kc.FourMomentumConstraints(p4, tlis, Vme);
373 kc.FourMomentumConstraints(p4, tlis, Vme);
374 m_kc.push_back(kc);
375 if((unsigned int) number != m_kc.size()-1)
376 std::cout << "wrong kinematic constraints index" << std::endl;
377}
378
379void KinematicFit::AddFourMomentum(int number, double etot ){
380
381 HepLorentzVector p4(0.0, 0.0, 0.0, etot);
382 std::vector<int> tlis;
383 tlis.clear();
384 KinematicConstraints kc;
385
386 for(int i = 0; i < numberWTrack(); i++) {
387 tlis.push_back(i);
388 }
389 HepSymMatrix Vme = HepSymMatrix (4,0);
390 Vme[3][3] = 2*m_espread*m_espread;
391 // kc.FourMomentumConstraints(p4, tlis, Vme);
392 kc.FourMomentumConstraints(p4, tlis, Vme);
393 m_kc.push_back(kc);
394 if((unsigned int) number != m_kc.size()-1)
395 std::cout << "wrong kinematic constraints index" << std::endl;
396}
397/*
398 void KinematicFit::AddPosition(int number, HepPoint3D xorigin, std::vector<int> tlis_V){
399 KinematicConstraints kc;
400 HepSymMatrix Vpe = HepSymMatrix(2,0);
401 HepSymMatrix Vclus = HepSymMatrix (3,0);
402 Vclus[0][0] = wTrackOrigin_V(tlis_V[0]).Ew()[4][4];
403 Vclus[0][1] = wTrackOrigin_V(tlis_V[0]).Ew()[4][5];
404 Vclus[1][1] = wTrackOrigin_V(tlis_V[0]).Ew()[5][5];
405 Vclus[2][2] = wTrackOrigin_V(tlis_V[0]).Ew()[6][6];
406 HepMatrix p(2,3,0);
407 p[0][0] = wTrackOrigin_V(tlis_V[0]).w()[1];
408 p[0][1] = -wTrackOrigin_V(tlis_V[0]).w()[0];
409 p[1][0] = wTrackOrigin_V(tlis_V[0]).w()[2];
410 p[1][2] = -wTrackOrigin_V(tlis_V[0]).w()[0];
411 Vpe = Vclus.similarity(p);
412
413 kc.PositionConstraints(xorigin, tlis_V, Vpe);
414 m_kc.push_back(kc);
415 if((unsigned int) number != m_kc.size()-1)
416 std::cout << "wrong kinematic constraints index" << std::endl;
417 }
418 */
419
420
421
422void KinematicFit::fit() {
423
424 KinematicConstraints kc;
425 int nc = m_nc;
426 int ntrk = numberWTrack();
427
428 TStopwatch timer;
429 timer.Start();
430
431 m_VD = HepSymMatrix(m_nc, 0);
432 m_VD = m_covOrigin.similarity(m_D);
433 timer.Stop();
434 m_cpu[1] += timer.CpuTime();
435 timer.Start();
436
437 int ifail;
438 m_VD.invert(ifail);
439
440 timer.Stop();
441 m_cpu[2] += timer.CpuTime();
442 timer.Start();
443
444 HepVector G(m_nc, 0);
445 G = m_D * (m_pOrigin - m_pInfit) + m_d;
446
447 m_KP = HepMatrix(NTRKPAR*m_nktrk, m_nc, 0);
448 m_KP = m_covOrigin * m_DT * m_VD;
449 m_chi = (m_VD.similarity(G.T()))[0][0];
450
451
452 timer.Stop();
453 m_cpu[3] += timer.CpuTime();
454 timer.Start();
455
456 m_pInfit = m_pOrigin - m_KP * G;
457
458 // ===== gamma dynamic adjustment====
459
460 timer.Stop();
461 m_cpu[4] += timer.CpuTime();
462
463 timer.Start();
464 if(m_dynamicerror ==1)
465 gda();
466
467 timer.Stop();
468 m_cpu[6] += timer.CpuTime();
469}
470
471
472
473
474
475void KinematicFit::upCovmtx() {
476
477 int ntrk = numberWTrack();
478 HepSymMatrix I(NTRKPAR*m_nktrk, 1);
479 m_covInfit = m_covOrigin.similarity(I - m_KP*m_D);
480 for(int i = 0; i < m_nktrk; i++) {
481 HepSymMatrix ew3 (3, 0);
482 HepSymMatrix ew6 (6, 0);
483 HepSymMatrix Ew1 (7, 0);
484 ew3 = m_covInfit.sub(i*NTRKPAR + 1, (i+1)*NTRKPAR);
485 for(int j = 0; j < NTRKPAR; j++) {
486 for(int k = 0; k < NTRKPAR; k++) {
487 ew6[j][k] = ew3[j][k];
488 }
489 }
490 for(int m = NTRKPAR; m < 6; m++) {
491 for(int n = NTRKPAR; n < 6; n++) {
492 ew6[m][n] = 0;//(wTrackOrigin(i).Ew())[m+1][n+1];
493 }
494 }
495 double px = p4Infit(i).px();
496 double py = p4Infit(i).py();
497 double pz = p4Infit(i).pz();
498 double m = p4Infit(i).m();
499 double e = p4Infit(i).e();
500
501 HepMatrix J(7, 6, 0);
502 J[0][0] = 1;
503 J[1][1] = 1;
504 J[2][2] = 1;
505 J[3][0] = px/e;
506 J[3][1] = py/e;
507 J[3][2] = pz/e;
508 J[4][3] = 1;
509 J[5][4] = 1;
510 J[6][5] = 1;
511 Ew1 = ew6.similarity(J);
512
513 HepVector W(7, 0);
514 for(int j=0; j<4; j++) {
515 W[j] = p4Infit(i)[j];
516 }
517 W[4] = wTrackOrigin(i).w()[4];
518 W[5] = wTrackOrigin(i).w()[5];
519 W[6] = wTrackOrigin(i).w()[6];
520
521
522 WTrackParameter wtrk = wTrackInfit(i);
523 wtrk.setEw(Ew1);
524 wtrk.setW(W);
525 setWTrackInfit(i, wtrk);
526 }
527
528
529}
530
531
532
533void KinematicFit::fit(int n) {
534
535 if(m_chisq.size() == 0) {
536 for(unsigned int i = 0; i < m_kc.size(); i++)
537 m_chisq.push_back(9999.0);
538 }
539 KinematicConstraints kc;
540 int nc = m_nc;
541 int ntrk = numberWTrack();
542
543 m_VD = HepSymMatrix(m_nc, 0);
544 m_VD = m_covOrigin.similarity(m_D);
545 int ifail;
546 m_VD.invert(ifail);
547 HepVector G(m_nc, 0);
548 G = m_D * (m_pOrigin - m_pInfit) + m_d;
549 m_KP = HepMatrix(NTRKPAR*m_nktrk, m_nc, 0);
550 m_KP = m_covOrigin * m_DT * m_VD;
551 m_chisq[n] = (m_VD.similarity(G.T()))[0][0];
552 m_pInfit = m_pOrigin - m_KP * G;
553
554}
555
556
557
558void KinematicFit::covMatrix(int n) {
559 KinematicConstraints kc = m_kc[n];
560 int nc = kc.nc();
561 int ntrk = (kc.Ltrk()).size();
562 HepSymMatrix Va0(7*ntrk, 0);
563 for (int i = 0; i < ntrk; i++) {
564 int itk = (kc.Ltrk())[i];
565 for(int j = 0; j < 7; j++)
566 for (int k = 0; k < 7; k++)
567 Va0[7*i+j][7*i+k] = (wTrackOrigin(itk).Ew())[j][k];
568
569 }
570 HepMatrix D(nc, 7*ntrk, 0);
571 int ncc = 0;
572 for(int j = 0; j < kc.nc(); j++) {
573 for(int l = 0; l < ntrk; l++) {
574 for(int k = 0; k < 7; k++)
575 D[ncc][7*l+k] = (kc.Dc()[l])[j][k];
576 }
577 ncc++;
578 }
579
580
581 HepSymMatrix VD(nc, 0);
582 VD = kc.VD()[0];
583 HepSymMatrix Va(7*ntrk, 0);
584 Va = Va0 - (VD.similarity(D.T())).similarity(Va0);
585 for(int i = 0; i < ntrk; i++) {
586 int itk = (kc.Ltrk())[i];
587 HepVector w(7, 0);
588 HepSymMatrix Ew(7, 0);
589 for(int j = 0; j < 7; j++) {
590 for(int k = 0; k < 7; k++)
591 Ew[j][k] = Va[7*i + j] [7*i + k];
592 }
593 WTrackParameter wtrk = wTrackInfit(itk);
594 wtrk.setEw(Ew);
595 setWTrackInfit(itk, wtrk);
596 }
597 m_kc[n] = kc;
598
599}
600
601
602
603
604bool KinematicFit::Fit() {
605 bool okfit = false;
606 TStopwatch timer;
607 m_nktrk = numberWTrack();
608 m_pOrigin = HepVector(m_nktrk * NTRKPAR, 0);
609 m_pInfit = HepVector(m_nktrk * NTRKPAR, 0);
610 m_covOrigin = HepSymMatrix(m_nktrk * NTRKPAR, 0);
611 m_covInfit = HepSymMatrix(m_nktrk * NTRKPAR, 0);
612 m_massvector = HepVector(m_nktrk, 0);
613 for(int i = 0; i < numberWTrack(); i++) {
614 setWTrackInfit(i, wTrackOrigin(i));
615 setPOrigin(i, (wTrackOrigin(i).w()).sub(1, NTRKPAR));
616 setPInfit(i, (wTrackOrigin(i).w()).sub(1, NTRKPAR));
617 setCovOrigin(i, (wTrackOrigin(i).Ew()).sub(1, NTRKPAR));
618 setMassvector(i, wTrackOrigin(i).mass());
619 }
620
621 //
622 // iteration
623 //
624 // cout<<"m_pInfit ="<<m_pInfit<<endl;
625 // cout<<"m_covOrigin="<<m_covOrigin<<endl;
626 // cout<<"m_massvector ="<<m_massvector<<endl;
627
628 std::vector<double> chisq;
629 chisq.clear();
630 int nc = 0;
631 for(int i = 0; i < m_kc.size(); i++)
632 nc += m_kc[i].nc();
633
634 m_D = HepMatrix(nc, m_nktrk * NTRKPAR, 0);
635 m_DT = HepMatrix(m_nktrk * NTRKPAR, nc, 0);
636 m_d = HepVector(nc, 0);
637
638 for(int it = 0; it < m_niter; it++) {
639
640 timer.Start();
641 m_nc = 0;
642 for(unsigned int i = 0; i < m_kc.size(); i++) {
643 KinematicConstraints kc = m_kc[i];
644 // std::vector<WTrackParameter> wlis;
645 // std::vector<WTrackParameter> wlis_V;
646 // wlis.clear();
647 // wlis_V.clear();
648 // for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
649 // int n = (kc.Ltrk())[j];
650 // WTrackParameter wtrk = wTrackInfit(n);
651 // if(m_espread!=0) wtrk = wTrackOrigin(n);
652 // wlis.push_back(wtrk);
653 // }
654 // for(unsigned int j = 0; j < (kc.Ltrk_V()).size(); j++) {
655 // int n = (kc.Ltrk_V())[j];
656 // WTrackParameter wtrk = wTrackInfit_V(n);
657 // wlis_V.push_back(wtrk);
658 // }
659 // kc.UpdateConstraints(wlis, wlis_V);
660 // m_kc[i] = kc;
661 //cout<<"wlis_V ="<<(wlis_V[0].w())[0]<<endl;
662 updateConstraints(kc);
663 // std::cout << "updata OK " << m_d << std::endl;
664 }
665 timer.Stop();
666 m_cpu[0] += timer.CpuTime();
667
668 fit();
669 chisq.push_back(m_chi);
670 if(it > 0) {
671 double delchi = chisq[it]- chisq[it-1];
672 if(fabs(delchi) < m_chiter)
673 break;
674 }
675 }
676 if(m_chi >= m_chicut) {
677
678 return okfit;
679 }
680 // update track parameter and its covariance matrix
681 // upTrkpar();
682 // covMatrix();
683 timer.Start();
684 //upCovmtx();
685 timer.Stop();
686 m_cpu[5] += timer.CpuTime();
687
688 okfit = true;
689
690 /*
691 for (int i = 0; i<numberWTrack(); i++){
692 if (wTrackOrigin(i).charge()==0) continue ;
693 HTrackParameter horigin = HTrackParameter(wTrackOrigin(i));
694 HTrackParameter hinfit = HTrackParameter(wTrackInfit(i));
695
696 HepVector a0 = horigin.hel();
697 HepVector a1 = hinfit.hel();
698 HepSymMatrix v0 = horigin.eHel();
699 HepSymMatrix v1 = hinfit.eHel();
700 HepVector pull(5,0);
701 for (int k=0; k<5; k++) {
702 pull[k] = (a0[k]-a1[k])/sqrt(abs(v0[k][k]-v1[k][k]));
703 }
704
705 WTrackParameter wtrk2 = wTrackInfit(i);
706 wtrk2.setPull(pull);
707 // for (int l=0;l<5; l++) {
708 //(wTrackInfit(i).pull())[l]=(wtrk2.pull())[l];
709 // }
710 setWTrackInfit(i, wtrk2);
711 }
712 */
713 /*/
714
715 for (int i = 0; i<numberWTrack_V(); i++){
716 //if (wTrackOrigin(i).charge()==0) continue ;
717 HTrackParameter horigin_V = HTrackParameter(wTrackOrigin_V(i));
718 HTrackParameter hinfit_V = HTrackParameter(wTrackInfit_V(i));
719
720 HepVector a0 = horigin.hel();
721 HepVector a1 = hinfit.hel();
722 HepSymMatrix v0 = horigin.eHel();
723 HepSymMatrix v1 = hinfit.eHel();
724 HepVector pull(5,0);
725 for (int k=0; k<5; k++) {
726 pull[k] = (a0[k]-a1[k])/sqrt(abs(v0[k][k]-v1[k][k]));
727 }
728
729 WTrackParameter wtrk2 = wTrackInfit(i);
730 wtrk2.setPull(pull);
731 // for (int l=0;l<5; l++) {
732 //(wTrackInfit(i).pull())[l]=(wtrk2.pull())[l];
733 // }
734 setWTrackInfit(i, wtrk2);
735 }
736 */
737
738 return okfit;
739}
740
741bool KinematicFit::Fit(int n) {
742 bool okfit = false;
743 if(n < 0 || (unsigned int)n >= m_kc.size()) return okfit;
744
745 m_nktrk = numberWTrack();
746 m_pOrigin = HepVector(m_nktrk * NTRKPAR, 0);
747 m_pInfit = HepVector(m_nktrk * NTRKPAR, 0);
748 m_covOrigin = HepSymMatrix(m_nktrk * NTRKPAR, 0);
749 m_covInfit = HepSymMatrix(m_nktrk * NTRKPAR, 0);
750 m_massvector = HepVector(m_nktrk * NTRKPAR, 0);
751 for(int i = 0; i < numberWTrack(); i++) {
752 setWTrackInfit(i, wTrackOrigin(i));
753 setPOrigin(i, (wTrackOrigin(i).w()).sub(1, NTRKPAR));
754 setPInfit(i, (wTrackOrigin(i).w()).sub(1, NTRKPAR));
755 setCovOrigin(i, (wTrackOrigin(i).Ew()).sub(1, NTRKPAR));
756 setMassvector(i, wTrackOrigin(i).mass());
757 }
758
759 //
760 // iteration loop
761 //
762
763 std::vector<double> chisq;
764 chisq.clear();
765
766 m_D = HepMatrix(m_kc[n].nc(), m_nktrk * NTRKPAR, 0);
767 m_DT = HepMatrix(m_nktrk * NTRKPAR, m_kc[n].nc(), 0);
768 m_d = HepVector(m_kc[n].nc(), 0);
769
770 for(int it = 0; it < m_niter; it++) {
771 m_nc = 0;
772 KinematicConstraints kc = m_kc[n];
773 updateConstraints(kc);
774 // m_kc[n] = kc;
775 fit(n);
776
777 chisq.push_back(m_chisq[n]);
778 if(it > 0) {
779 double delchi = chisq[it]- chisq[it-1];
780 if(fabs(delchi) < m_chiter)
781 break;
782 }
783 }
784
785 if(m_chisq[n] >= m_chicut) return okfit;
786 // ====update cov====
787 // upCovmtx();
788 okfit = true;
789 return okfit;
790}
791
792
793void KinematicFit::BuildVirtualParticle( int n) {
794//
795// q = p1 + p2 + ... + pn
796//
797 upCovmtx();
798 KinematicConstraints kc = m_kc[n];
799 int ntrk = (kc.Ltrk()).size();
800 int charge = 0;
801 HepVector w(7, 0);
802 HepSymMatrix ew(7, 0);
803 HepMatrix dwdp(7, 7, 0);
804 dwdp[0][0] = 1;
805 dwdp[1][1] = 1;
806 dwdp[2][2] = 1;
807 dwdp[3][3] = 1;
808 dwdp[4][4] = 1;
809 dwdp[5][5] = 1;
810 dwdp[6][6] = 1;
811 for (int i = 0; i < ntrk; i++) {
812 int itk = (kc.Ltrk())[i];
813 charge += wTrackInfit(itk).charge();
814 w[0] = w[0] + wTrackInfit(itk).w()[0];
815 w[1] = w[1] + wTrackInfit(itk).w()[1];
816 w[2] = w[2] + wTrackInfit(itk).w()[2];
817 w[3] = w[3] + wTrackInfit(itk).w()[3];
818 w[4] = 0.0;//
819 w[5] = 0.0;// set virtual particle's vertex at (0,0,0)
820 w[6] = 0.0;//
821 ew = ew + (wTrackInfit(itk).Ew()).similarity(dwdp); // the vertex matrix of this particles is not correct, because we do not use vertex information in kinematicfit, so ...
822 }
823 double m = sqrt(w[3]*w[3] - w[0]*w[0] - w[1]*w[1] - w[2]*w[2]);
824 WTrackParameter vwtrk;
825 vwtrk.setCharge(charge);
826 vwtrk.setW(w);
827 vwtrk.setEw(ew);
828 vwtrk.setMass(m);
829 m_virtual_wtrk.push_back(vwtrk);
830}
831
832
833
834void KinematicFit::gda(){
835 for(int i = 0; i < numberWTrack(); i++) {
836
837 if ((wTrackOrigin(i)).type() == 2 ){
838 int n ;
839 for(int j = 0; j < numberGammaShape(); j++) {
840 if(i==GammaShapeList(j)) n = j;
841 }
842 HepSymMatrix Ew(NTRKPAR, 0);
843 HepLorentzVector p1 = p4Infit(i);
844 double eorigin = pOrigin(i)[3];
845 //cout<<"p1 ="<<p1<<endl;
846 HepMatrix dwda1(NTRKPAR, 3, 0);
847 dwda1[0][0] = -p1.py();
848 dwda1[1][0] = p1.px();
849 dwda1[0][1] = p1.px()*p1.pz()/p1.perp();
850 dwda1[1][1] = p1.py()*p1.pz()/p1.perp();
851 dwda1[2][1] = -p1.perp();
852 dwda1[0][2] = p1.px()/p1.rho();
853 dwda1[1][2] = p1.py()/p1.rho();
854 dwda1[2][2] = p1.pz()/p1.rho();
855 dwda1[3][2] = p1.rho()/p1.e();
856 HepSymMatrix emcmea1(3, 0);
857 double pk = p1[3];
858 emcmea1[0][0] = GammaShapeValue(n).getdphi() * GammaShapeValue(n).getdphi();
859 emcmea1[1][1] = GammaShapeValue(n).getdthe() * GammaShapeValue(n).getdthe();
860 emcmea1[2][2] = GammaShapeValue(n).de(eorigin,pk) * GammaShapeValue(n).de(eorigin,pk);
861 Ew = emcmea1.similarity(dwda1);
862 //cout<<"emcmea1 ="<<emcmea1<<endl;
863 //cout<<"Ew ="<<Ew<<endl;
864 setCovOrigin(i, Ew);
865 }
866 }
867}
868
869
870HepVector KinematicFit::pull(int n){
871 upCovmtx();
872
873 if (wTrackOrigin(n).charge()!=0){
874 HepVector W(6,0);
875 HepSymMatrix Ew(6,0);
876 HepVector W1(7,0);
877 HepSymMatrix Ew1(7,0);
878 WTrackParameter wtrk = wTrackOrigin(n);
879 // W = wTrackOrigin(n).w();
880 // Ew = wTrackOrigin(n).Ew();
881 //cout<<"===Origin status==="<<endl;
882 //cout<<"W = "<<W<<endl;
883 //cout<<"Ew ="<<Ew<<endl;
884 for(int i=0; i<3; i++) {
885 W[i] = pInfit(n)[i];
886 }
887 W[3] = wTrackOrigin(n).w()[4];
888 W[4] = wTrackOrigin(n).w()[5];
889 W[5] = wTrackOrigin(n).w()[6];
890 for(int j=0; j<3; j++) {
891 for(int k=0; k<3; k++) {
892 Ew[j][k] = covInfit(n)[j][k];
893 }
894 }
895
896 for(int j=3; j<6; j++) {
897 for(int k=3; k<6; k++) {
898 Ew[j][k] = wTrackOrigin(n).Ew()[j+1][k+1];
899 }
900 }
901 //
902 // define J matrix to transfer 3 parameters to 4 parameters
903 //
904 double px = p4Infit(n).px();
905 double py = p4Infit(n).py();
906 double pz = p4Infit(n).pz();
907 double e = p4Infit(n).e();
908 HepMatrix J(7, 6, 0);
909 J[0][0] = 1;
910 J[1][1] = 1;
911 J[2][2] = 1;
912 J[3][0] = px/e;
913 J[3][1] = py/e;
914 J[3][2] = pz/e;
915 J[4][3] = 1;
916 J[5][4] = 1;
917 J[6][5] = 1;
918 W1 = J * W;
919 Ew1 = Ew.similarity(J) ;
920
921
922
923 // cout<<"===Infiddt status==="<<endl;
924 // cout<<"p4 ="<<p4Infit(n)<<endl;
925 // cout<<"W ="<<wTrackOrigin(n).w()<<endl;
926 // cout<<"W1 ="<<W1<<endl;
927 // cout<<"Ew ="<<wTrackOrigin(n).Ew()<<endl;
928 // cout<<"Ew1 ="<<Ew1<<endl;
929
930 wtrk.setW(W1);
931 wtrk.setEw(Ew1);
932 setWTrackInfit(n, wtrk);
933 HTrackParameter horigin = HTrackParameter(wTrackOrigin(n));
934 HTrackParameter hinfit = HTrackParameter(wTrackInfit(n));
935
936 HepVector a0 = horigin.hel();
937 HepVector a1 = hinfit.hel();
938 HepSymMatrix v0 = horigin.eHel();
939 HepSymMatrix v1 = hinfit.eHel();
940 HepVector pull(11,0);
941 for (int k=0; k<5; k++) {
942 pull[k] = (a0[k]-a1[k])/sqrt(abs(v0[k][k]-v1[k][k]));
943 // cout<<"pull ["<<k<<"] ="<<pull[k]<<endl;
944 }
945 for (int l=5; l<9; l++) {
946 pull[l] = (wTrackOrigin(n).w()[l-5] - wTrackInfit(n).w()[l-5])/sqrt(abs(wTrackOrigin(n).Ew()[l-5][l-5] - wTrackInfit(n).Ew()[l-5][l-5]));
947 // cout<<"pull ["<<l<<"] ="<<pull[l]<<endl;
948 }
949
950 // pull[9] = wTrackOrigin(n).w()[3] - wTrackInfit(n).w()[3];
951 // pull[10] =1/sqrt(abs(wTrackOrigin(n).Ew()[3][3] - wTrackInfit(n).Ew()[3][3]));
952 return pull;
953 }else {
954 HepVector pull(3,0);
955 for (int m=0; m<3; m++) {
956 pull[m] = (wTrackOrigin(n).w()[m] - wTrackInfit(n).w()[m])/sqrt(abs(wTrackOrigin(n).Ew()[m][m] - wTrackInfit(n).Ew()[m][m]));
957 }
958 return pull;
959 }
960}
961
962
963void KinematicFit::updateConstraints(KinematicConstraints k ) {
964 KinematicConstraints kc = k;
965
966 // std::vector<HepLorentzVector> wlis;
967 // std::vector<WTrackParameter> wlis_V;
968 // wlis.clear();
969 // wlis_V.clear();
970 // for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
971 // int n = (kc.Ltrk())[j];
972 // HepVec wtrk = wTrackInfit(n);
973 // if(m_espread!=0) wtrk = wTrackOrigin(n);
974
975 // wlis.push_back(p4Infit(n));
976 // }
977 // for(unsigned int j = 0; j < (kc.Ltrk_V()).size(); j++) {
978 // int n = (kc.Ltrk_V())[j];
979 // WTrackParameter wtrk = wTrackInfit_V(n);
980 // wlis_V.push_back(wtrk);
981 // }
982
983
984 int type = kc.Type();
985 switch(type) {
986 case Resonance: {
987 //
988 // E^2 - px^2 - py^2 - pz^2 = mres^2
989 //
990 double mres = kc.mres();
991 HepLorentzVector pmis;
992 for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
993 int n = (kc.Ltrk())[j];
994
995 pmis = pmis + p4Infit(n);
996 }
997 // for(unsigned int i = 0; i < wlis_V.size(); i++)
998 // pmis = pmis + wlis_V[i].p();
999 for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
1000 int n = (kc.Ltrk())[j];
1001 HepMatrix Dc(1, NTRKPAR, 0);
1002 Dc[0][0] = -2 * pmis.px() + 2 * pmis.e() * p4Infit(n).px() / p4Infit(n).e();
1003 Dc[0][1] = -2 * pmis.py() + 2 * pmis.e() * p4Infit(n).py() / p4Infit(n).e();
1004 Dc[0][2] = -2 * pmis.pz() + 2 * pmis.e() * p4Infit(n).pz() / p4Infit(n).e();
1005 // Dc[0][3] = 2 * pmis.e();
1006 setD(m_nc,n,Dc);
1007 setDT(n, m_nc, Dc.T());
1008 }
1009
1010 // for(unsigned int i = 0; i < wlis_V.size(); i++) {
1011 // HepMatrix Dc_V(1, 7, 0);
1012 // Dc_V[0][0] = -2 * pmis.px();
1013 // Dc_V[0][1] = -2 * pmis.py();
1014 // Dc_V[0][2] = -2 * pmis.pz();
1015 // Dc_V[0][3] = 2 * pmis.e();
1016 // m_Dc_V.push_back(Dc_V);
1017 // }
1018
1019
1020 HepVector dc(1, 0);
1021 dc[0] = pmis.m2() - mres * mres;
1022 m_d[m_nc] = dc[0];
1023 m_nc+=1;
1024 // std::cout << "etot = " << dc[0] <<" , " << mres<< std::endl;
1025 // m_dc.push_back(dc);
1026 // HepVector lambda(1, 0);
1027 // m_lambda.push_back(lambda);
1028 // HepSymMatrix vd(1, 0);
1029 // m_VD.push_back(vd);
1030 // HepSymMatrix V6 = m_Vre;
1031 // m_Vm.push_back(V6);
1032
1033 break;
1034 }
1035 case TotalEnergy: {
1036 //
1037 // E - Etot = 0
1038 //
1039 double etot = kc.etot();
1040 HepLorentzVector pmis;
1041 for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++){
1042 int n = (kc.Ltrk())[j];
1043 pmis = pmis + p4Infit(n);
1044 }
1045 // for(unsigned int i = 0; i < wlis.size(); i++)
1046 // pmis = pmis + wlis[i].p();
1047 for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
1048 int n = (kc.Ltrk())[j];
1049 HepMatrix Dc(1, NTRKPAR, 0);
1050 Dc[0][0] = p4Infit(n).px()/p4Infit(n).e();
1051 Dc[0][1] = p4Infit(n).py()/p4Infit(n).e();
1052 Dc[0][2] = p4Infit(n).pz()/p4Infit(n).e();
1053 // Dc[0][3] = 1.0;
1054 setD(m_nc,n,Dc);
1055 setDT(n, m_nc, Dc.T());
1056 }
1057 HepVector dc(1, 0);
1058 dc[0] = pmis.e() - etot;
1059 m_d[m_nc] = dc[0];
1060 m_nc+=1;
1061 // m_dc.push_back(dc);
1062 // HepVector lambda(1, 0);
1063 // m_lambda.push_back(lambda);
1064 // HepSymMatrix vd(1, 0);
1065 // m_VD.push_back(vd);
1066 break;
1067 }
1068 case TotalMomentum: {
1069 //
1070 // sqrt(px^2+py^2+pz^2) - ptot = 0
1071 //
1072 double ptot = kc.ptot();
1073 HepLorentzVector pmis;
1074 for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
1075 int n = (kc.Ltrk())[j];
1076 pmis = pmis + p4Infit(n);
1077 }
1078
1079 // for(unsigned int i = 0; i < wlis.size(); i++)
1080 // pmis = pmis + wlis[i].p();
1081 for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
1082 int n = (kc.Ltrk())[j];
1083 HepMatrix Dc(1, NTRKPAR, 0);
1084 Dc[0][0] = pmis.px()/pmis.rho();
1085 Dc[0][1] = pmis.py()/pmis.rho();
1086 Dc[0][2] = pmis.pz()/pmis.rho();
1087 setD(m_nc,n,Dc);
1088 setDT(n, m_nc, Dc.T());
1089 // m_Dc.push_back(Dc);
1090 }
1091 HepVector dc(1, 0);
1092 dc[0] = pmis.rho() - ptot;
1093 m_d[m_nc] = dc[0];
1094 m_nc+=1;
1095 // m_dc.push_back(dc);
1096 // HepVector lambda(1, 0);
1097 // m_lambda.push_back(lambda);
1098 // HepSymMatrix vd(1, 0);
1099 // m_VD.push_back(vd);
1100 break;
1101 }
1102 /*
1103 case kc.typePoint(): {
1104 HepPoint3D point = kc.point();
1105 double flightpath;
1106 HepMatrix p(2,3,0);
1107
1108 HepSymMatrix Vp(2,0);
1109 for(unsigned int i = 0; i < wlis.size(); i++) {
1110 HepMatrix Dc(2, 7, 0);
1111 m_Dc.push_back(Dc);
1112 }
1113 // wlis_V.size() should be 1
1114 for(unsigned int i = 0; i < wlis_V.size(); i++) {
1115 // HepMatrix Dc(3, 7, 0);
1116 // m_Dc.push_back(Dc);
1117 HepMatrix Dc_V(2, 7, 0);
1118 HepSymMatrix Vclus(3, 0);
1119 for(unsigned int j = 0; j < 3; j++){
1120 for(unsigned int k = j; k < 3; k++){
1121 Vclus[j][k] = wlis_V[i].Ew()[j+4][k+4];
1122 }
1123 }
1124 cout<<"Vclus ="<<Vclus<<endl;
1125 p[0][0] = wlis_V[i].w()[1];
1126 p[0][1] = -wlis_V[i].w()[0];
1127 p[1][0] = wlis_V[i].w()[2];
1128 p[1][2] = -wlis_V[i].w()[0];
1129 Vp = Vclus.similarity(p);
1130 cout<<"Vp ="<<Vp<<endl;
1131 Dc_V[0][0] = -(wlis_V[i].w()[5] - point[1]);
1132 Dc_V[0][1] = wlis_V[i].w()[4] - point[0];
1133 Dc_V[0][4] = wlis_V[i].w()[1];
1134 Dc_V[0][5] = -wlis_V[i].w()[0];
1135
1136 Dc_V[1][0] = -(wlis_V[i].w()[6] - point[2]);
1137 Dc_V[1][2] = wlis_V[i].w()[4] - point[0];
1138 Dc_V[1][4] = wlis_V[i].w()[2];
1139 Dc_V[1][6] = -wlis_V[i].w()[0];
1140
1141 // HepMatrix q_x(3,1,0) ;
1142 // q_x[0][0] = 1;
1143 // HepMatrix deltaX_x(3,1,0) ;
1144 // deltaX_x[0][0] = 1;
1145 // HepMatrix p1_x = -(q_x.T()*Vclus.inverse(iver)*deltaX)*p2 + p1*(q_x.T()*Vclus.inverse(iver)*q + q.T()*Vclus.inverse(iver)*q_x);
1146 // HepMatrix p2_x = p2*p2;
1147
1148 // HepMatrix q_y(3,1,0) ;
1149 // q_y[1][0] = 1;
1150 // HepMatrix deltaX_y(3,1,0);
1151 // deltaX_y[1][0] = 1;
1152 // HepMatrix p1_y = -(q_y.T()*Vclus.inverse(iver)*deltaX)*p2 + p1*(q_y.T()*Vclus.inverse(iver)*q + q.T()*Vclus.inverse(iver)*q_y);
1153 // HepMatrix p2_y = p2*p2;
1154
1155 // HepMatrix q_z(3,1,0);
1156 // q_z[2][0] = 1;
1157 // HepMatrix deltaX_z(3,1,0);
1158 // deltaX_z[2][0] = 1;
1159 // Hw()[0]<<endl;
1160 //cout<<"dc[0] ="<<dc[0]<<endl;
1161 dc[1] = (wlis_V[i].w()[4] - point[0])*wlis_V[i].w()[2] - (wlis_V[i].w()[6] - point[2])*wlis_V[i].w()[0];
1162 m_dc.push_back(dc);
1163 }
1164 HepSymMatrix V3 = Vp;
1165 m_Vm.push_back(V3);
1166 break;
1167 }
1168 */
1169
1170 case ThreeMomentum: {
1171 //
1172 // px - p3x = 0
1173 // py - p3y = 0
1174 // pz - p3z = 0
1175 //
1176 Hep3Vector p3 = kc.p3();
1177 HepLorentzVector pmis;
1178 for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
1179 int n = (kc.Ltrk())[j];
1180 pmis = pmis + p4Infit(n);
1181 }
1182 // for(unsigned int i = 0; i < wlis.size(); i++)
1183 // pmis = pmis + wlis[i].p();
1184 for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
1185 int n = (kc.Ltrk())[j];
1186
1187 HepMatrix Dc(kc.nc(), NTRKPAR, 0);
1188 Dc[0][0] = 1.0;
1189 Dc[1][1] = 1.0;
1190 Dc[2][2] = 1.0;
1191 HepVector dc(kc.nc(), 0);
1192 dc[0] = pmis.px() - p3.x();
1193 dc[1] = pmis.py() - p3.y();
1194 dc[2] = pmis.pz() - p3.z();
1195 for(int i = 0; i < kc.nc(); i++) {
1196 setD(m_nc+i, n, Dc.sub(i+1, i+1, 1, NTRKPAR));
1197 setDT(n, m_nc+i, (Dc.sub(i+1, i+1, 1, NTRKPAR)).T());
1198 m_d[m_nc+i] = dc[i];
1199 }
1200 // m_Dc.push_back(Dc);
1201 }
1202 m_nc += 3;
1203
1204 // HepVector dc(3, 0);
1205 // dc[0] = pmis.px() - p3.x();
1206 // dc[1] = pmis.py() - p3.y();
1207 // dc[2] = pmis.pz() - p3.z();
1208 // m_dc.push_back(dc);
1209 // HepVector lambda(3, 0);
1210 // m_lambda.push_back(lambda);
1211 // HepSymMatrix vd(3, 0);
1212 // m_VD.push_back(vd);
1213 break;
1214 }
1215 case EqualMass: {
1216 //
1217 // (E_1 ^2 - Px_1 ^2 - Py_1 ^2 - Pz_1 ^2) - (E_2 ^2 - Px_2 ^2 - Py_2 ^2 - Pz_2 ^2) = 0
1218 //
1219
1220 int isiz = (kc.numEqual())[0];
1221 HepLorentzVector pmis1, pmis2;
1222 for(int n = 0; n < isiz; n++) {
1223 int n1 = (kc.Ltrk())[n];
1224 pmis1 = pmis1 + p4Infit(n1);
1225 }
1226 int jsiz = (kc.numEqual())[1];
1227 for(int n = 0; n < jsiz; n++){
1228 int n2 = (kc.Ltrk())[n+isiz];
1229 pmis2 = pmis2 + p4Infit(n2);
1230 }
1231 for(int i = 0; i < isiz; i++) {
1232 int n1 = (kc.Ltrk())[i];
1233 HepMatrix Dc(1, NTRKPAR, 0);
1234 Dc[0][0] = -2 * pmis1.px() + 2 * pmis1.e() * p4Infit(n1).px() / p4Infit(n1).e();
1235 Dc[0][1] = -2 * pmis1.py() + 2 * pmis1.e() * p4Infit(n1).py() / p4Infit(n1).e();
1236 Dc[0][2] = -2 * pmis1.pz() + 2 * pmis1.e() * p4Infit(n1).pz() / p4Infit(n1).e();
1237 // Dc[0][3] = 2 * pmis1.e();
1238 setD(m_nc,n1,Dc);
1239 setDT(n1,m_nc,Dc.T());
1240 }
1241 for(int i = 0; i < jsiz; i++) {
1242 int n2 = (kc.Ltrk())[i+isiz];
1243 HepMatrix Dc(1, NTRKPAR, 0);
1244 Dc[0][0] = 2 * pmis2.px() - 2 * pmis2.e() * p4Infit(n2).px() / p4Infit(n2).e();
1245 Dc[0][1] = 2 * pmis2.py() - 2 * pmis2.e() * p4Infit(n2).py() / p4Infit(n2).e();
1246 Dc[0][2] = 2 * pmis2.pz() - 2 * pmis2.e() * p4Infit(n2).pz() / p4Infit(n2).e();
1247 Dc[0][3] = -2 * pmis2.e();
1248 setD(m_nc,n2,Dc);
1249 setDT(n2,m_nc,Dc.T());
1250 }
1251 // int isiz_V = m_nequal_V[0];
1252 // HepLorentzVector pmis1_V, pmis2_V;
1253 // if(isiz_V > 0){
1254 // for(int n = 0; n < isiz_V; n++) {
1255 // pmis1_V = pmis1_V + wlis_V[n].p();
1256 // }
1257 // }
1258 // int jsiz_V = m_nequal_V[1];
1259 // if(jsiz_V > 0) {
1260 // for(int n = 0; n < jsiz_V; n++)
1261 // pmis2_V = pmis2_V + wlis_V[isiz_V+n].p();
1262 // }
1263
1264 // for(int i = 0; i < isiz_V; i++) {
1265 // HepMatrix Dc_V(1, 7, 0);
1266 // Dc_V[0][0] = -2 * pmis1_V.px();
1267 // Dc_V[0][1] = -2 * pmis1_V.py();
1268 // Dc_V[0][2] = -2 * pmis1_V.pz();
1269 // Dc_V[0][3] = 2 * pmis1_V.e();
1270 // m_Dc_V.push_back(Dc_V);
1271 // }
1272 // for(int i = isiz_V; i < isiz_V+jsiz_V; i++) {
1273 // HepMatrix Dc_V(1, 7, 0);
1274 // Dc_V[0][0] = 2 * pmis2_V.px();
1275 // Dc_V[0][1] = 2 * pmis2_V.py();
1276 // Dc_V[0][2] = 2 * pmis2_V.pz();
1277 // Dc_V[0][3] = -2 * pmis2_V.e();
1278 // m_Dc_V.push_back(Dc_V);
1279 // }
1280 HepVector dc(1, 0);
1281 dc[0] = pmis1.m2() - pmis2.m2();// + pmis1_V.m2() - pmis2_V.m2();
1282 m_d[m_nc] = dc[0];
1283
1284 m_nc+=1;
1285 // m_dc.push_back(dc);
1286 // HepVector lambda(1, 0);
1287 // m_lambda.push_back(lambda);
1288 // HepSymMatrix vd(1, 0);
1289 // m_VD.push_back(vd);
1290 // HepSymMatrix V2 = m_Vne;
1291 // m_Vm.push_back(V2);
1292 // std::cout <<"m_Vm[0] ="<<m_Vm[0]<<std::endl;
1293 break;
1294 }
1295 case FourMomentum:
1296 default: {
1297 //
1298 // px - p4x = 0
1299 // py - p4y = 0
1300 // pz - p4z = 0
1301 // e - p4e = 0
1302 //
1303 HepLorentzVector p4 = kc.p4();
1304 HepLorentzVector pmis;
1305 for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
1306 int n = (kc.Ltrk())[j];
1307 pmis = pmis + p4Infit(n);
1308 }
1309 // HepLorentzVector pmis_V;
1310 for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
1311 int n = (kc.Ltrk())[j];
1312 HepMatrix Dc(kc.nc(), NTRKPAR, 0);
1313 Dc[0][0] = 1.0;
1314 Dc[1][1] = 1.0;
1315 Dc[2][2] = 1.0;
1316 Dc[3][0] = p4Infit(n).px() / p4Infit(n).e();
1317 Dc[3][1] = p4Infit(n).py() / p4Infit(n).e();
1318 Dc[3][2] = p4Infit(n).pz() / p4Infit(n).e();
1319 // Dc[3][3] = 1.0;
1320
1321 // m_Dc.push_back(Dc);
1322 HepVector dc(kc.nc(), 0);
1323 dc[0] = pmis.px() - p4.px();
1324 dc[1] = pmis.py() - p4.py();
1325 dc[2] = pmis.pz() - p4.pz();
1326 dc[3] = pmis.e() - p4.e();
1327 for(int i = 0; i < kc.nc(); i++) {
1328 setD(m_nc+i, n, Dc.sub(i+1, i+1, 1, NTRKPAR));
1329 setDT(n, m_nc+i, (Dc.sub(i+1, i+1, 1, NTRKPAR)).T());
1330 m_d[m_nc+i] = dc[i];
1331 }
1332 }
1333 m_nc += 4;
1334 // for(unsigned int i = 0; i < wlis_V.size(); i++)
1335 // pmis_V = pmis_V + wlis_V[i].p();
1336 // for(unsigned int i = 0; i < wlis_V.size(); i++) {
1337 // HepMatrix Dc_V(4, 7, 0);
1338 // Dc_V[0][0] = 1.0;
1339 // Dc_V[1][1] = 1.0;
1340 // Dc_V[2][2] = 1.0;
1341 // Dc_V[3][3] = 1.0;
1342 // m_Dc_V.push_back(Dc_V);
1343 // }
1344
1345 // HepVector dc(4, 0);
1346 // dc[0] = pmis.px() + pmis_V.px() - p4.px();
1347 // dc[1] = pmis.py() + pmis_V.py() - p4.py();
1348 // dc[2] = pmis.pz() + pmis_V.pz() - p4.pz();
1349 // dc[3] = pmis.e() + pmis_V.e() - p4.e();
1350 // m_dc.push_back(dc);
1351
1352 // HepSymMatrix V1 = m_Vme;
1353 // m_Vm.push_back(V1);
1354 // HepVector lambda(4, 0);
1355 // m_lambda.push_back(lambda);
1356 // HepSymMatrix vd(4, 0);
1357 // m_VD.push_back(vd);
1358
1359 break;
1360 }
1361 }
1362}
1363
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
mass
double mass
Definition: CosmicGenerator.cxx:138
n
const Int_t n
Definition: DataBase/tau_mode.c:65
etot
Double_t etot
Definition: DataBase/tau_mode.c:7
I
const DifComplex I
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
dc
int dc[18]
Definition: EvtPycont.cc:41
HepPoint3D
HepGeom::Point3D< double > HepPoint3D
Definition: Gam4pikp.cxx:37
n2
int n2
Definition: SD0Tag.cxx:55
n1
int n1
Definition: SD0Tag.cxx:54
HTrackParameter::hel
HepVector hel() const
Definition: HTrackParameter.h:49
HTrackParameter::eHel
HepSymMatrix eHel() const
Definition: HTrackParameter.h:51
KinematicConstraints::VD
std::vector< HepSymMatrix > VD()
Definition: KinematicConstraints.h:44
KinematicConstraints::FourMomentumConstraints
void FourMomentumConstraints(const HepLorentzVector p4, std::vector< int > tlis, HepSymMatrix Vme)
Definition: KinematicConstraints.cxx:74
KinematicConstraints::Ltrk
std::vector< int > Ltrk()
Definition: KinematicConstraints.h:45
KinematicConstraints::ResonanceConstraints
void ResonanceConstraints(const double mres, std::vector< int > tlis, HepSymMatrix Vre)
Definition: KinematicConstraints.cxx:34
KinematicConstraints::ThreeMomentumConstraints
void ThreeMomentumConstraints(const Hep3Vector p3, std::vector< int > tlis)
Definition: KinematicConstraints.cxx:64
KinematicConstraints::p3
Hep3Vector p3() const
Definition: KinematicConstraints.h:65
KinematicConstraints::TotalEnergyConstraints
void TotalEnergyConstraints(const double etot, std::vector< int > tlis)
Definition: KinematicConstraints.cxx:44
KinematicConstraints::numEqual
std::vector< int > numEqual()
Definition: KinematicConstraints.h:74
KinematicConstraints::p4
HepLorentzVector p4() const
Definition: KinematicConstraints.h:66
KinematicConstraints::Dc
std::vector< HepMatrix > Dc()
Definition: KinematicConstraints.h:40
KinematicConstraints::TotalMomentumConstraints
void TotalMomentumConstraints(const double ptot, std::vector< int > tlis)
Definition: KinematicConstraints.cxx:54
KinematicConstraints::EqualMassConstraints
void EqualMassConstraints(std::vector< int > tlis1, std::vector< int > tlis2, HepSymMatrix Vne)
Definition: KinematicConstraints.cxx:85
KinematicFit::instance
static KinematicFit * instance()
Definition: KinematicFit.cxx:20
KinematicFit::BuildVirtualParticle
void BuildVirtualParticle(int number)
Definition: KinematicFit.cxx:793
KinematicFit::AddResonance
void AddResonance(int number, double mres, std::vector< int > tlis)
Definition: KinematicFit.cxx:143
KinematicFit::AddEqualMass
void AddEqualMass(int number, std::vector< int > tlis1, std::vector< int > tlis2)
Definition: KinematicFit.cxx:322
KinematicFit::chisq
double chisq() const
Definition: KinematicFit.h:150
KinematicFit::AddThreeMomentum
void AddThreeMomentum(int number, Hep3Vector p3)
Definition: KinematicFit.cxx:335
KinematicFit::AddFourMomentum
void AddFourMomentum(int number, HepLorentzVector p4)
Definition: KinematicFit.cxx:354
KinematicFit::pull
HepVector pull(int n)
Definition: KinematicFit.cxx:870
KinematicFit::AddTotalEnergy
void AddTotalEnergy(int number, double etot, std::vector< int > lis)
Definition: KinematicFit.cxx:311
KinematicFit::AddTotalMomentum
void AddTotalMomentum(int number, double ptot, std::vector< int > lis)
Definition: KinematicFit.cxx:228
TrackPool::AddList
std::vector< int > AddList(int n1)
Definition: TrackPool.cxx:483
TrackPool::wTrackInfit
std::vector< WTrackParameter > wTrackInfit() const
Definition: TrackPool.h:73
TrackPool::setVBeamPosition
void setVBeamPosition(const HepSymMatrix VBeamPosition)
Definition: TrackPool.h:144
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::numberGammaShape
int numberGammaShape() const
Definition: TrackPool.h:99
TrackPool::clearWTrackList
void clearWTrackList()
Definition: TrackPool.h:113
TrackPool::clearMappositionA
void clearMappositionA()
Definition: TrackPool.h:125
TrackPool::GammaShapeValue
std::vector< GammaShape > GammaShapeValue() const
Definition: TrackPool.h:95
TrackPool::wTrackOrigin
std::vector< WTrackParameter > wTrackOrigin() const
Definition: TrackPool.h:72
TrackPool::GammaShapeList
std::vector< int > GammaShapeList() const
Definition: TrackPool.h:96
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::setBeamPosition
void setBeamPosition(const HepPoint3D BeamPosition)
Definition: TrackPool.h:143
TrackPool::setWTrackInfit
void setWTrackInfit(const int n, const WTrackParameter wtrk)
Definition: TrackPool.h:106
WTrackParameter::setEw
void setEw(const HepSymMatrix &Ew)
Definition: WTrackParameter.h:53
WTrackParameter::setMass
void setMass(const double mass)
Definition: WTrackParameter.h:55
WTrackParameter::setCharge
void setCharge(const int charge)
Definition: WTrackParameter.h:54
WTrackParameter::setW
void setW(const HepVector &w)
Definition: WTrackParameter.h:51
Alignment::NTRKPAR
const int NTRKPAR
Definition: Alignment.h:49